/* nag_zgges (f08xnc) Example Program. * * Copyright 2011 Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #include #include #include int main(void) { /* Scalars */ Complex alph, bet, z; double norma, normb, normd, norme, eps, small; Integer i, j, n, sdim, pda, pdb, pdc, pdd, pde, pdvsl, pdvsr; Integer exit_status = 0; /* Arrays */ Complex *a = 0, *alpha = 0, *b = 0, *beta = 0, *c = 0; Complex *d = 0, *e = 0, *vsl = 0, *vsr = 0; Nag_LeftVecsType jobvsl; Nag_RightVecsType jobvsr; char nag_enum_arg[40]; /* Nag Types */ NagError fail; Nag_OrderType order; #ifdef NAG_COLUMN_MAJOR #define A(I, J) a[(J-1)*pda + I - 1] #define B(I, J) b[(J-1)*pdb + I - 1] order = Nag_ColMajor; #else #define A(I, J) a[(I-1)*pda + J - 1] #define B(I, J) b[(I-1)*pdb + J - 1] order = Nag_RowMajor; #endif INIT_FAIL(fail); printf("nag_zgges (f08xnc) Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n]"); scanf("%ld%*[^\n]", &n); if (n < 0) { printf("Invalid n\n"); exit_status = 1; goto END; } scanf(" %s%*[^\n]", nag_enum_arg); /* nag_enum_name_to_value(x04nac). * Converts NAG enum member name to value */ jobvsl = (Nag_LeftVecsType) nag_enum_name_to_value(nag_enum_arg); scanf(" %s%*[^\n]", nag_enum_arg); jobvsr = (Nag_RightVecsType) nag_enum_name_to_value(nag_enum_arg); pdvsl = (jobvsl==Nag_LeftVecs?n:1); pdvsr = (jobvsr==Nag_RightVecs?n:1); pda = n; pdb = n; pdc = n; pdd = pda; pde = pdb; /* Allocate memory */ if (!(a = NAG_ALLOC(n*n, Complex)) || !(b = NAG_ALLOC(n*n, Complex)) || !(c = NAG_ALLOC(n*n, Complex)) || !(d = NAG_ALLOC(n*n, Complex)) || !(e = NAG_ALLOC(n*n, Complex)) || !(alpha = NAG_ALLOC(n, Complex)) || !(beta = NAG_ALLOC(n, Complex)) || !(vsl = NAG_ALLOC(pdvsl*pdvsl, Complex)) || !(vsr = NAG_ALLOC(pdvsr*pdvsr, Complex))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Read in the matrices A and B */ for (i = 1; i <= n; ++i) for (j = 1; j <= n; ++j) scanf(" ( %lf , %lf )", &A(i, j).re, &A(i, j).im); scanf("%*[^\n]"); for (i = 1; i <= n; ++i) for (j = 1; j <= n; ++j) scanf(" ( %lf , %lf )", &B(i, j).re, &B(i, j).im); scanf("%*[^\n]"); /* nag_zge_norm (f16uac): Find norms of input matrices A and B. */ nag_zge_norm(order, Nag_OneNorm, n, n, a, pda, &norma, &fail); nag_zge_norm(order, Nag_OneNorm, n, n, b, pdb, &normb, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dge_norm (f16rac).\n%s\n", fail.message); exit_status = 1; goto END; } /* Copy Matrices A and B to D and E respectively: nag_zge_copy (f16tfc), * Complex valued general matrix copy. */ nag_zge_copy(order, Nag_NoTrans, n, n, a, pda, d, pdd, &fail); nag_zge_copy(order, Nag_NoTrans, n, n, b, pdb, e, pde, &fail); /* nag_gen_complx_mat_print_comp (x04dbc): Print matrix A */ fflush(stdout); nag_gen_complx_mat_print_comp(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, n, a, pda, Nag_BracketForm, "%6.2f", "Matrix A", Nag_IntegerLabels, 0, Nag_IntegerLabels, 0, 80, 0, 0, &fail); printf("\n"); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_complx_mat_print_comp (x04dbc).\n%s\n", fail.message); exit_status = 1; goto END; } /* nag_gen_complx_mat_print_comp (x04dbc): Print matrix B */ fflush(stdout); nag_gen_complx_mat_print_comp(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, n, b, pdb, Nag_BracketForm, "%6.2f", "Matrix B", Nag_IntegerLabels, 0, Nag_IntegerLabels, 0, 80, 0, 0, &fail); printf("\n"); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_complx_mat_print_comp (x04dbc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Find the generalized Schur form using nag_zgges (f08xnc). */ nag_zgges(order, jobvsl, jobvsr, Nag_NoSortEigVals, Nag_FALSE, n, a, pda, b, pdb, &sdim, alpha, beta, vsl, pdvsl, vsr, pdvsr, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_zgges (f08xnc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Print Eigenvalues */ /* nag_real_safe_small_number (x02amc), nag_machine_precision (x02ajc) */ eps = nag_machine_precision; small = nag_real_safe_small_number; printf("\n Eigenvalues\n"); for (j = 0; j < n; ++j) { printf("%2ld ", j+1); if (nag_complex_abs(alpha[j]) * small >= nag_complex_abs(beta[j])) { printf(" numerically infinite or undetermined, alpha = " "(%9.4f, %9.4f), beta = (%9.4f, %9.4f)\n", alpha[j].re, alpha[j].im, beta[j].re, beta[j].im); } else { z = nag_complex_divide(alpha[j], beta[j]); printf(" (%11.4e, %11.4e)\n", z.re, z.im); } } /* Check generalized Schur Form by reconstruction of Schur vectors are * available. */ if (jobvsl==Nag_NotLeftVecs || jobvsr==Nag_NotRightVecs) { /* Cannot check factorization by reconstruction Schur vectors. */ goto END; } /* Reconstruct A as Q*S*Z^H and subtract from original (D) using the steps * C = Q (Q in vsl) using nag_zge_copy (f16tfc). * C = C*S (S in a, upper triangular) using nag_ztrmm (f16zfc). * D = D - C*Z^H (Z in vsr) using nag_zgemm (f16zac). */ nag_zge_copy(order, Nag_NoTrans, n, n, vsl, pdvsl, c, pdc, &fail); alph = nag_complex(1.0,0.0); /* nag_ztrmm (f16zfc) Triangular complex matrix-matrix multiply. */ nag_ztrmm(order, Nag_RightSide, Nag_Upper, Nag_NoTrans, Nag_NonUnitDiag, n, n, alph, a, pda, c, pdc, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ztrmm (f16zfc).\n%s\n", fail.message); exit_status = 1; goto END; } alph = nag_complex(-1.0,0.0); bet = nag_complex(1.0,0.0); nag_zgemm(order, Nag_NoTrans, Nag_ConjTrans, n, n, n, alph, c, pdc, vsr, pdvsr, bet, d, pdd, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_zgemm (f16zac).\n%s\n", fail.message); exit_status = 1; goto END; } /* Reconstruct B as Q*T*Z^H and subtract from original (D) using the steps * Q = Q*T (Q in vsl, T in b, upper triangular) using nag_ztrmm (f16zfc). * E = E - Q*Z^H (Z in vsr) using nag_zgemm (f16zac). */ alph = nag_complex(1.0,0.0); /* nag_ztrmm (f16zfc) Triangular complex matrix-matrix multiply. */ nag_ztrmm(order, Nag_RightSide, Nag_Upper, Nag_NoTrans, Nag_NonUnitDiag, n, n, alph, b, pdb, vsl, pdvsl, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ztrmm (f16zfc).\n%s\n", fail.message); exit_status = 1; goto END; } alph = nag_complex(-1.0,0.0); bet = nag_complex(1.0,0.0); nag_zgemm(order, Nag_NoTrans, Nag_ConjTrans, n, n, n, alph, vsl, pdvsl, vsr, pdvsr, bet, e, pde, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_zgemm (f16zac).\n%s\n", fail.message); exit_status = 1; goto END; } /* Compute norm of difference matrices D and E using * nag_zge_norm (f16uac). */ nag_zge_norm(order, Nag_OneNorm, n, n, d, pdd, &normd, &fail); nag_zge_norm(order, Nag_OneNorm, n, n, e, pde, &norme, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_zge_norm (f16uac).\n%s\n", fail.message); exit_status = 1; goto END; } if (MAX(normd,norme) > pow(eps,0.8)*MAX(norma,normb)) { printf("The norm of the error in the reconstructed matrices is greater " "than expected.\nThe Schur factorization has failed.\n"); exit_status = 1; } END: if (a) NAG_FREE(a); if (b) NAG_FREE(b); if (c) NAG_FREE(c); if (d) NAG_FREE(d); if (e) NAG_FREE(e); if (alpha) NAG_FREE(alpha); if (beta) NAG_FREE(beta); if (vsl) NAG_FREE(vsl); if (vsr) NAG_FREE(vsr); return exit_status; }