f08wbce.c

/* nag_lapackeig_dggevx (f08wbc) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<math.h>
#include<nag.h>
#include<stdio.h>
intmain(void){
/* Scalars */
Complexeig,eigl,eigr;
doubleabnorm,abnrm,bbnrm,eps,sign,small,tol;
Integeri,ihi,ilo,j,k,n,pda,pdb,pdvl,pdvr;
Integerverbose=0;
Integerexit_status=0;
/* Arrays */
double*a=0,*alphai=0,*alphar=0,*b=0,*beta=0;
double*lscale=0,*rconde=0,*rcondv=0,*rscale=0;
double*vl=0,*vr=0;
charnag_enum_arg[40];
/* Nag Types */
NagErrorfail;
Nag_OrderTypeorder;
Nag_LeftVecsTypejobvl;
Nag_RightVecsTypejobvr;
Nag_RCondTypesense;
#ifdef NAG_COLUMN_MAJOR
#define A(I, J) a[(J - 1) * pda + I - 1]
#define B(I, J) b[(J - 1) * pdb + I - 1]
#define VL(I, J) vl[(J - 1) * pdvl + I - 1]
#define VR(I, J) vr[(J - 1) * pdvr + 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]
#define VL(I, J) vl[(I - 1) * pdvl + J - 1]
#define VR(I, J) vr[(I - 1) * pdvr + J - 1]
order=Nag_RowMajor;
#endif
INIT_FAIL(fail);
printf("nag_lapackeig_dggevx (f08wbc) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%"NAG_IFMT"%*[^\n]",&n);
if(n<0){
printf("Invalid n\n");
exit_status=1;
gotoEND;
}
scanf(" %39s%*[^\n]",nag_enum_arg);
/* nag_enum_name_to_value (x04nac).
 * Converts NAG enum member name to value
 */
jobvl=(Nag_LeftVecsType)nag_enum_name_to_value(nag_enum_arg);
scanf(" %39s%*[^\n]",nag_enum_arg);
jobvr=(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);
pda=n;
pdb=n;
pdvl=(jobvl==Nag_LeftVecs?n:1);
pdvr=(jobvr==Nag_RightVecs?n:1);
/* Allocate memory */
if(!(a=NAG_ALLOC(n*n,double))||!(alphai=NAG_ALLOC(n,double))||
!(alphar=NAG_ALLOC(n,double))||!(b=NAG_ALLOC(n*n,double))||
!(beta=NAG_ALLOC(n,double))||!(lscale=NAG_ALLOC(n,double))||
!(rconde=NAG_ALLOC(n,double))||!(rcondv=NAG_ALLOC(n,double))||
!(rscale=NAG_ALLOC(n,double))||
!(vl=NAG_ALLOC(pdvl*pdvl,double))||
!(vr=NAG_ALLOC(pdvr*pdvr,double))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
/* 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]");
/* Solve the generalized eigenvalue problem using nag_lapackeig_dggevx
 * (f08wbc). */
nag_lapackeig_dggevx(order,Nag_BalanceBoth,jobvl,jobvr,sense,n,a,pda,
b,pdb,alphar,alphai,beta,vl,pdvl,vr,pdvr,&ilo,
&ihi,lscale,rscale,&abnrm,&bbnrm,rconde,rcondv,
&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dggevx (f08wbc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_machine_real_safe (x02amc), nag_machine_precision (x02ajc) */
eps=nag_machine_precision;
small=nag_machine_real_safe;
if(abnrm==0.0)
abnorm=ABS(bbnrm);
elseif(bbnrm==0.0)
abnorm=ABS(abnrm);
elseif(ABS(abnrm)>=ABS(bbnrm))
abnorm=ABS(abnrm)*sqrt(1.0+(bbnrm/abnrm)*(bbnrm/abnrm));
else
abnorm=ABS(bbnrm)*sqrt(1.0+(abnrm/bbnrm)*(abnrm/bbnrm));
tol=eps*abnorm;
/* Print out eigenvalues and vectors and associated condition
 * number and bounds.
 */
for(j=0;j<n;++j){
/* Print out information on the j-th eigenvalue */
printf("\n");
if((fabs(alphar[j])+fabs(alphai[j]))*small>=fabs(beta[j])){
printf("Eigenvalue %2"NAG_IFMT" is numerically infinite or "
"undetermined\n",
j+1);
printf("alpha = (%13.4e, %13.4e), beta = %13.4e\n",alphar[j],alphai[j],
beta[j]);
}elseif(alphai[j]==0.0){
printf("Eigenvalue %2"NAG_IFMT" = %13.4e\n",j+1,
alphar[j]/beta[j]);
}else{
eig.re=alphar[j]/beta[j],eig.im=alphai[j]/beta[j];
printf("Eigenvalue %2"NAG_IFMT" = (%13.4e, %13.4e)\n",j+1,eig.re,
eig.im);
}
if(verbose){
if(sense==Nag_RCondEigVals||sense==Nag_RCondBoth){
printf("\n Reciprocal condition number = %10.1e\n",rconde[j]);
if(rconde[j]>0.0)
printf(" Error bound = %10.1e\n",tol/rconde[j]);
else
printf(" Error bound is infinite\n");
}
}
printf("\n\n");
/* Normalize and print out information on the j-th eigenvector(s) */
if(jobvl==Nag_LeftVecs)
printf("%21s%8s","Left Eigenvector","");
if(jobvr==Nag_RightVecs)
printf("%21s","Right Eigenvector");
printf(" %2"NAG_IFMT"\n",j+1);
if(alphai[j]==0.0)
for(i=1;i<=n;++i){
if(jobvl==Nag_LeftVecs)
printf("%7s%13.4e%12s","",VL(i,j+1)/VL(n,j+1),"");
if(jobvr==Nag_RightVecs)
printf("%7s%13.4e","",VR(i,j+1)/VR(n,j+1));
printf("\n");
}
else{
k=(alphai[j]>0.0?j+1:j);
sign=(alphai[j]>0.0?1.0:-1.0);
if(jobvl==Nag_LeftVecs)
eigl=nag_complex_create(VL(n,k),VL(n,k+1));
if(jobvr==Nag_RightVecs)
eigr=nag_complex_create(VR(n,k),VR(n,k+1));
for(i=1;i<=n;++i){
if(jobvl==Nag_LeftVecs){
eig=nag_complex_divide(nag_complex_create(VL(i,k),VL(i,k+1)),
eigl);
printf(" (%13.4e,%13.4e) ",eig.re,sign*eig.im);
}
if(jobvr==Nag_RightVecs){
eig=nag_complex_divide(nag_complex_create(VR(i,k),VR(i,k+1)),
eigr);
printf(" (%13.4e,%13.4e)",eig.re,sign*eig.im);
}
printf("\n");
}
}
if(verbose){
if(sense==Nag_RCondEigVecs||sense==Nag_RCondBoth){
printf("\n Reciprocal condition number = %10.1e\n",rcondv[j]);
if(rcondv[j]>0.0)
printf(" Error bound = %10.1e\n\n",tol/rcondv[j]);
else
printf(" Error bound is infinite\n\n");
}
}
}
END:
NAG_FREE(a);
NAG_FREE(alphai);
NAG_FREE(alphar);
NAG_FREE(b);
NAG_FREE(beta);
NAG_FREE(lscale);
NAG_FREE(rconde);
NAG_FREE(rcondv);
NAG_FREE(rscale);
NAG_FREE(vl);
NAG_FREE(vr);
returnexit_status;
}