/* nag_dggesx (f08xbc) Example Program. * * Copyright 2011 Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #include #include static Nag_Boolean NAG_CALL delctg(const double ar, const double ai, const double b); int main(void) { /* Scalars */ double dg_a, dg_b, eps, norma, normb, normd, norme, small, 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_RCondType sense; #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 sort = Nag_SortEigVals; sense = Nag_RCondBoth; 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; 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); scanf(" %s%*[^\n]", nag_enum_arg); sense = (Nag_RCondType) nag_enum_name_to_value(nag_enum_arg); pda = n; pdb = n; pdc = n; pdd = pda; pde = pdb; pdvsl = (jobvsl==Nag_LeftVecs?n:1); pdvsr = (jobvsr==Nag_RightVecs?n:1); /* 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_FrobeniusNorm, n, n, a, pda, &norma, &fail); 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; } norma = sqrt(norma*norma + normb*normb); /* 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); 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 Matrix A. */ 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 (x04cac): Print Matrix B. */ 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_dggesx (f08xbc). */ nag_dggesx(order, jobvsl, jobvsr, sort, delctg, sense, n, a, pda, b, pdb, &sdim, alphar, alphai, beta, vsl, pdvsl, vsr, pdvsr, rconde, rcondv, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dggesx (f08xbc).\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; tol = eps * norma; /* 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]); } if (sense==Nag_RCondEigVals || sense==Nag_RCondBoth) { /* Print out the reciprocal condition number and error bound */ printf("\n"); printf("For the selected eigenvalues,\nthe reciprocals of projection " "norms onto the deflating subspaces are\n"); printf(" for left subspace, rcond = %10.1e\n for right subspace, rcond = " "%10.1e\n\n", rconde[0], rconde[1]); printf(" asymptotic error bound = %10.1e\n\n", tol / rconde[0]); } /* 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); 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; } normd = sqrt(normd*normd + norme*norme); if (normd > pow(eps,0.8)*norma) { 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 (sense==Nag_RCondEigVecs || sense==Nag_RCondBoth) { /* Print out the reciprocal condition numbers and error bound. */ printf("For the left and right deflating subspaces,\n"); printf("reciprocal condition numbers are:\n"); printf(" for left subspace, rcond = %10.1e\n for right subspace, rcond = " "%10.1e\n\n", rcondv[0], rcondv[1]); printf(" approximate error bound = %10.1e\n", tol / rcondv[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 (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; } Nag_Boolean NAG_CALL delctg(const double ar, const double ai, const double b) { /* Boolean function delctg for use with nag_dggesx (f08xbc) * Returns the value Nag_TRUE if the imaginary part of the eigenvalue * (ar + ai*i)/B is zero, i.e. the eigenvalue is real. */ return (ai == 0.0?Nag_TRUE: Nag_FALSE); }