/* nag_dgges (f08xac) Example Program. * * Copyright 2011 Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static Nag_Boolean NAG_CALL delctg(const double ar, const double ai, const double b); #ifdef __cplusplus } #endif int main(void) { /* Scalars */ Integer exit_status = 0; Integer i, j, n, sdim, pda, pdb, pdc, pdd, pde, pdvsl, pdvsr; double dg_a, dg_b, eps, small, norma, normb, normd, norme; /* Arrays */ double *a = 0, *alphai = 0, *alphar = 0, *b = 0, *beta = 0; double *c = 0, *d = 0, *e = 0, *vsl = 0, *vsr = 0; char nag_enum_arg[40]; /* Nag Types */ NagError fail; Nag_OrderType order; Nag_LeftVecsType jobvsl; Nag_RightVecsType jobvsr; #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_dgges (f08xac) 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, double)) || !(b = NAG_ALLOC(n*n, double)) || !(c = NAG_ALLOC(n*n, double)) || !(d = NAG_ALLOC(n*n, double)) || !(e = NAG_ALLOC(n*n, double)) || !(alphai = NAG_ALLOC(n, double)) || !(alphar = NAG_ALLOC(n, double)) || !(beta = NAG_ALLOC(n, double)) || !(vsl = NAG_ALLOC(pdvsl*pdvsl, double)) || !(vsr = NAG_ALLOC(pdvsr*pdvsr, double))) { 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", &A(i, j)); scanf("%*[^\n]"); for (i = 1; i <= n; ++i) for (j = 1; j <= n; ++j) scanf("%lf", &B(i, j)); scanf("%*[^\n]"); /* nag_dge_norm (f16rac): Find norms of input matrices A and B. */ nag_dge_norm(order, Nag_OneNorm, n, n, a, pda, &norma, &fail); nag_dge_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 matrices D and E using nag_dge_copy (f16qfc), * real valued general matrix copy. * The copies will be used as comparison against reconstructed matrices. */ nag_dge_copy(order, Nag_NoTrans, n, n, a, pda, d, pdd, &fail); if (fail.code == NE_NOERROR) { nag_dge_copy(order, Nag_NoTrans, n, n, b, pdb, e, pde, &fail); } if (fail.code != NE_NOERROR) { printf("Error from nag_dge_copy (f16qfc).\n%s\n", fail.message); exit_status = 1; goto END; } /* nag_gen_real_mat_print (x04cac): Print Matrices A and B. */ fflush(stdout); nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, n, a, pda, "Matrix A", 0, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n"); fflush(stdout); nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, n, b, pdb, "Matrix B", 0, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\n"); fflush(stdout); /* Find the generalized Schur form using nag_dgges (f08xac). */ nag_dgges(order, jobvsl, jobvsr, Nag_SortEigVals, delctg, n, a, pda, b, pdb, &sdim, alphar, alphai, beta, vsl, pdvsl, vsr, pdvsr, &fail); if (fail.code != NE_NOERROR && fail.code != NE_SCHUR_REORDER_SELECT) { printf("Error from nag_dgges (f08xac).\n%s\n", fail.message); exit_status = 1; goto END; } printf("Number of sorted eigenvalues = %4ld\n\n", sdim); /* nag_real_safe_small_number (x02amc), nag_machine_precision (x02ajc) */ eps = nag_machine_precision; small = nag_real_safe_small_number; /* Print the eigenvalues */ printf(" Eigenvalues\n"); for (j = 0; j < n; ++j) { printf("%2ld ", j+1); if ((fabs(alphar[j]) + fabs(alphai[j])) * small >= fabs(beta[j])) printf(" infinite or undetermined, alpha = (%11.4e, %11.4e), " "beta = %11.4e\n", alphar[j], alphai[j], beta[j]); else if (alphai[j] == 0.0) printf(" %12.4e\n", alphar[j]/beta[j]); else printf(" (%11.4e, %11.4e)\n", alphar[j]/beta[j], alphai[j]/beta[j]); } /* 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^T and subtract from original (D) using the steps * C = Q*S (Q in vsl, S in a) using nag_dgemm (f16yac). * Note: not nag_dtrmm since S may not be strictly triangular. * D = D - C*Z^T (Z in vsr) using nag_dgemm (f16yac). */ dg_a = 1.0; dg_b = 0.0; nag_dgemm(order, Nag_NoTrans, Nag_NoTrans, n, n, n, dg_a, vsl, pdvsl, a, pda, dg_b, c, pdc, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dgemm (f16yac).\n%s\n", fail.message); exit_status = 1; goto END; } dg_a = -1.0; dg_b = 1.0; nag_dgemm(order, Nag_NoTrans, Nag_Trans, n, n, n, dg_a, c, pdc, vsr, pdvsr, dg_b, d, pdd, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dgemm (f16yac).\n%s\n", fail.message); exit_status = 1; goto END; } /* Reconstruct B as Q*T*Z^T and subtract from original (E) using the steps * C = Q*T (Q in vsl, T in b) using nag_dgemm (f16yac). * E = E - C*Z^T (Z in vsr) using nag_dgemm (f16yac). */ dg_a = 1.0; dg_b = 0.0; nag_dgemm(order, Nag_NoTrans, Nag_NoTrans, n, n, n, dg_a, vsl, pdvsl, b, pdb, dg_b, c, pdc, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dgemm (f16yac).\n%s\n", fail.message); exit_status = 1; goto END; } dg_a = -1.0; dg_b = 1.0; nag_dgemm(order, Nag_NoTrans, Nag_Trans, n, n, n, dg_a, c, pdc, vsr, pdvsr, dg_b, e, pde, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dgemm (f16yac).\n%s\n", fail.message); exit_status = 1; goto END; } /* nag_dge_norm (f16rac): Find norms of difference matrices D and E. */ nag_dge_norm(order, Nag_OneNorm, n, n, d, pdd, &normd, &fail); nag_dge_norm(order, Nag_OneNorm, n, n, e, pde, &norme, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dge_norm (f16rac).\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; goto END; } if (fail.code == NE_SCHUR_REORDER_SELECT) printf("*** Note that rounding errors mean that leading eigenvalues in the" " generalized\n Schur form no longer satisfy delctg = Nag_TRUE\n\n"); 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 (alphai) NAG_FREE(alphai); if (alphar) NAG_FREE(alphar); if (beta) NAG_FREE(beta); if (vsl) NAG_FREE(vsl); if (vsr) NAG_FREE(vsr); return exit_status; } #undef B #undef A static Nag_Boolean NAG_CALL delctg(const double ar, const double ai, const double b) { /* Logical function delctg for use with nag_dgges (f08xac). * * Returns the value Nag_TRUE if the eigenvalue is real. */ return (ai == 0.0?Nag_TRUE:Nag_FALSE); }