f08pbce.c

/* nag_lapackeig_dgeesx (f08pbc) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<math.h>
#include<nag.h>
#include<stdio.h>
#ifdef __cplusplus
extern"C"{
#endif
staticNag_BooleanNAG_CALLselect_fun(constdoublewr,constdoublewi);
#ifdef __cplusplus
}
#endif
intmain(void){
/* Scalars */
doublealpha,anorm,beta,eps,norm,rconde,rcondv;
Integeri,j,n,pda,pdc,pdd,pdvs,sdim;
Integerexit_status=0;
/* Arrays */
double*a=0,*c=0,*d=0,*vs=0,*wi=0,*wr=0;
/* Nag Types */
NagErrorfail;
Nag_OrderTypeorder;
#ifdef NAG_COLUMN_MAJOR
#define A(I, J) a[(J - 1) * pda + I - 1]
order=Nag_ColMajor;
#else
#define A(I, J) a[(I - 1) * pda + J - 1]
order=Nag_RowMajor;
#endif
INIT_FAIL(fail);
printf("nag_lapackeig_dgeesx (f08pbc) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%"NAG_IFMT"%*[^\n]",&n);
if(n<0){
printf("Invalid n\n");
exit_status=1;
returnexit_status;
}
pda=n;
pdc=n;
pdd=n;
pdvs=n;
/* Allocate memory */
if(!(a=NAG_ALLOC(n*n,double))||!(c=NAG_ALLOC(n*n,double))||
!(d=NAG_ALLOC(n*n,double))||!(vs=NAG_ALLOC(n*n,double))||
!(wi=NAG_ALLOC(n,double))||!(wr=NAG_ALLOC(n,double))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
/* Read in the matrix A */
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
scanf("%lf",&A(i,j));
scanf("%*[^\n]");
/* Copy A to D: nag_blast_dge_copy (f16qfc),
 * real valued general matrix copy.
 */
nag_blast_dge_copy(order,Nag_NoTrans,n,n,a,pda,d,pdd,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_blast_dge_copy (f16qfc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_blast_dge_norm (f16rac): Find norm of matrix A for use later
 * in relative error test.
 */
nag_blast_dge_norm(order,Nag_OneNorm,n,n,a,pda,&anorm,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_blast_dge_norm (f16rac).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_file_print_matrix_real_gen (x04cac): Print Matrix A. */
fflush(stdout);
nag_file_print_matrix_real_gen(order,Nag_GeneralMatrix,Nag_NonUnitDiag,n,
n,a,pda,"Matrix A",0,&fail);
printf("\n");
if(fail.code!=NE_NOERROR){
printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n",
fail.message);
exit_status=1;
gotoEND;
}
/* Find the Schur factorization of A using nag_lapackeig_dgeesx (f08pbc). */
nag_lapackeig_dgeesx(order,Nag_Schur,Nag_SortEigVals,select_fun,
Nag_RCondBoth,n,a,pda,&sdim,wr,wi,vs,pdvs,
&rconde,&rcondv,&fail);
if(fail.code!=NE_NOERROR&&fail.code!=NE_SCHUR_REORDER_SELECT){
printf("Error from nag_lapackeig_dgeesx (f08pbc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Reconstruct A from Schur Factorization Z*T*Trans(Z) where T is upper
 * triangular and stored in A. This can be done using the following steps:
 * i. C = Z*T (nag_blast_dgemm, f16yac),
 * ii. D = D-C*trans(Z) (nag_blast_dgemm, f16yac).
 */
alpha=1.0;
beta=0.0;
nag_blast_dgemm(order,Nag_NoTrans,Nag_NoTrans,n,n,n,alpha,vs,pdvs,a,
pda,beta,c,pdc,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_blast_dgemm (f16yac).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_blast_dgemm (f16yac):
 * Compute D = A - C*Z^T.
 */
alpha=-1.0;
beta=1.0;
nag_blast_dgemm(order,Nag_NoTrans,Nag_Trans,n,n,n,alpha,c,pdc,vs,
pdvs,beta,d,pdd,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_blast_dgemm (f16yac).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_blast_dge_norm (f16rac): Find norm of difference matrix D and print
 * warning if it is too large relative to norm of A.
 */
nag_blast_dge_norm(order,Nag_OneNorm,n,n,d,pdd,&norm,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_blast_dge_norm (f16rac).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Get the machine precision, using nag_machine_precision (x02ajc) */
eps=nag_machine_precision;
if(norm>pow(eps,0.8)*MAX(anorm,1.0)){
printf("||A-(Z*T*Z^T)||/||A|| is larger than expected.\n"
"Schur factorization has failed.\n");
exit_status=1;
gotoEND;
}
/* Print details on eigenvalues */
printf("Number of eigenvalues for which select is true = %4"NAG_IFMT"\n\n",
sdim);
if(fail.code==NE_SCHUR_REORDER_SELECT){
printf(" ** Note that rounding errors mean that leading eigenvalues in the"
" Schur form\n no longer satisfy select(lambda) = Nag_TRUE\n\n");
}else{
printf("The selected eigenvalues are:\n");
for(i=0;i<sdim;i++)
printf("%3"NAG_IFMT" (%13.4e, %13.4e)\n",i+1,wr[i],wi[i]);
}
/* Print out the reciprocal condition numbers */
printf("\nReciprocal of projection norm onto the invariant subspace\n");
printf("%26sfor the selected eigenvalues rconde = %8.1e\n\n","",rconde);
printf("Reciprocal condition number for the invariant subspace rcondv = "
"%8.1e\n\n",
rcondv);
/* Compute the approximate asymptotic error bound on the average absolute
 * error of the selected eigenvalues given by eps*norm(A)/rconde.
 */
printf("Approximate asymptotic error bound for selected eigenvalues = "
"%8.1e\n\n",
eps*anorm/rconde);
/* Compute an approximate asymptotic bound on the maximum angular error in
 * the computed invariant subspace given by eps*norm(A)/rcondv
 */
printf("Approximate asymptotic error bound for the invariant subspace = "
"%8.1e\n",
eps*anorm/rcondv);
END:
NAG_FREE(a);
NAG_FREE(c);
NAG_FREE(d);
NAG_FREE(vs);
NAG_FREE(wi);
NAG_FREE(wr);
returnexit_status;
}
staticNag_BooleanNAG_CALLselect_fun(constdoublear,constdoubleai){
/* Boolean function select for use with nag_lapackeig_dgees (f08pac)
 * Returns the value Nag_TRUE if the eigenvalue is real and positive
 */
return(ar>0.0&&ai==0.0?Nag_TRUE:Nag_FALSE);
}