/* nag_dggesx (f08xbc) Example Program. * * Copyright 2013 Numerical Algorithms Group. * * Mark 24, 2013. */ #include #include #include #include #include #include #include #include #ifdef __cplusplus extern "C" { #endif static Nag_Boolean NAG_CALL selctg(const double ar, const double ai, const double b); #ifdef __cplusplus } #endif int main(void) { /* Scalars */ double abnorm, dg_a, dg_b, eps, norma, normb, normd, norme, tol; Integer i, j, n, sdim, pda, pdb, pdc, pdd, pde, pdvsl, pdvsr; Integer exit_status = 0; /* Arrays */ double *a = 0, *alphai = 0, *alphar = 0, *b = 0, *beta = 0; double *c = 0, *d = 0, *e =0, *vsl = 0, *vsr = 0; double rconde[2], rcondv[2]; char nag_enum_arg[40]; /* Nag Types */ NagError fail; Nag_OrderType order; Nag_LeftVecsType jobvsl; Nag_RightVecsType jobvsr; Nag_SortEigValsType sort = Nag_SortEigVals; Nag_RCondType sense = Nag_RCondBoth; #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_dggesx (f08xbc) 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; return exit_status; } scanf(" %39s%*[^\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(" %39s%*[^\n]", nag_enum_arg); jobvsr = (Nag_RightVecsType) nag_enum_name_to_value(nag_enum_arg); scanf(" %39s%*[^\n]", nag_enum_arg); sense = (Nag_RCondType) 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 = n; pde = n; /* 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]"); /* 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) { printf("Error from nag_dge_copy (f16qfc).\n%s\n", fail.message); exit_status = 1; goto END; } 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_dge_norm (f16rac): Find norms of input matrices A and B. */ nag_dge_norm(order, Nag_FrobeniusNorm, n, n, a, pda, &norma, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dge_norm (f16rac).\n%s\n", fail.message); exit_status = 1; goto END; } nag_dge_norm(order, Nag_FrobeniusNorm, 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; } /* 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); printf("\n"); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message); exit_status = 1; goto END; } fflush(stdout); nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, n, n, b, pdb, "Matrix B", 0, &fail); printf("\n"); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message); exit_status = 1; goto END; } /* Find the generalized Schur form using nag_dggesx (f08xbc). */ nag_dggesx(order, jobvsl, jobvsr, sort, selctg, sense, n, a, pda, b, pdb, &sdim, alphar, alphai, beta, vsl, pdvsl, vsr, pdvsr, rconde, rcondv, &fail); if (fail.code != NE_NOERROR && fail.code != NE_SCHUR_REORDER_SELECT) { printf("Error from nag_dggesx (f08xbc).\n%s\n", fail.message); exit_status = 1; goto END; } /* 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_FrobeniusNorm, n, n, d, pdd, &normd, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dge_norm (f16rac).\n%s\n", fail.message); exit_status = 1; goto END; } nag_dge_norm(order, Nag_FrobeniusNorm, 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; } /* Get the machine precision, using nag_machine_precision (x02ajc) */ eps = nag_machine_precision; 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; } /* Print details on eigenvalues */ printf("Number of sorted eigenvalues = %4ld\n\n", sdim); 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 selctg = Nag_TRUE" "\n\n"); } else { printf("The selected eigenvalues are:\n"); for (i=0;i 0.0 && ai == 0.0 && b != 0.0 ? Nag_TRUE : Nag_FALSE); }