f08ubce.c

/* nag_lapackeig_dsbgvx (f08ubc) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<nag.h>
#include<stdio.h>
intmain(void){
/* Scalars */
doubleabstol,vl,vu;
Integerexit_status=0,il=1,iu=1;
Integeri,j,ka,kb,m,n,pdab,pdbb,pdq,pdz,zsize;
/* Arrays */
double*ab=0,*bb=0,*q=0,*w=0,*z=0;
Integer*index=0;
charnag_enum_arg[40];
/* Nag Types */
NagErrorfail;
Nag_OrderTypeorder;
Nag_UploTypeuplo;
Nag_JobTypejob;
#ifdef NAG_COLUMN_MAJOR
#define AB_UPPER(I, J) ab[(J - 1) * pdab + ka + I - J]
#define AB_LOWER(I, J) ab[(J - 1) * pdab + I - J]
#define BB_UPPER(I, J) bb[(J - 1) * pdbb + kb + I - J]
#define BB_LOWER(I, J) bb[(J - 1) * pdbb + I - J]
order=Nag_ColMajor;
#else
#define AB_UPPER(I, J) ab[(I - 1) * pdab + J - I]
#define AB_LOWER(I, J) ab[(I - 1) * pdab + ka + J - I]
#define BB_UPPER(I, J) bb[(I - 1) * pdbb + J - I]
#define BB_LOWER(I, J) bb[(I - 1) * pdbb + kb + J - I]
order=Nag_RowMajor;
#endif
INIT_FAIL(fail);
printf("nag_lapackeig_dsbgvx (f08ubc) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%"NAG_IFMT"%"NAG_IFMT"%"NAG_IFMT"%*[^\n]",&n,&ka,&kb);
if(n<0||ka<kb||kb<0){
printf("Invalid n, ka or kb\n");
exit_status=1;
gotoEND;
}
scanf(" %39s%*[^\n]",nag_enum_arg);
/* nag_enum_name_to_value (x04nac).
 * Converts NAG enum member name to value
 */
uplo=(Nag_UploType)nag_enum_name_to_value(nag_enum_arg);
scanf(" %39s%*[^\n]",nag_enum_arg);
job=(Nag_JobType)nag_enum_name_to_value(nag_enum_arg);
if(job==Nag_EigVals){
zsize=1;
pdz=1;
}else{
zsize=n*n;
pdz=n;
}
pdab=ka+1;
pdbb=kb+1;
pdq=n;
/* Allocate memory */
if(!(ab=NAG_ALLOC((ka+1)*n,double))||
!(bb=NAG_ALLOC((kb+1)*n,double))||
!(q=NAG_ALLOC(n*n,double))||!(w=NAG_ALLOC(n,double))||
!(z=NAG_ALLOC(zsize,double))||!(index=NAG_ALLOC(n,Integer))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
/* Read the lower and upper bounds of the interval to be searched. */
scanf("%lf%lf%*[^\n]",&vl,&vu);
/* Read the triangular parts of the matrices A and B from data file */
if(uplo==Nag_Upper){
for(i=1;i<=n;++i)
for(j=i;j<=MIN(i+ka,n);++j)
scanf("%lf",&AB_UPPER(i,j));
scanf("%*[^\n]");
for(i=1;i<=n;++i)
for(j=i;j<=MIN(i+kb,n);++j)
scanf("%lf",&BB_UPPER(i,j));
}else{
for(i=1;i<=n;++i)
for(j=MAX(1,i-ka);j<=i;++j)
scanf("%lf",&AB_LOWER(i,j));
scanf("%*[^\n]");
for(i=1;i<=n;++i)
for(j=MAX(1,i-kb);j<=i;++j)
scanf("%lf",&BB_LOWER(i,j));
}
scanf("%*[^\n]");
/* Use the default absolute error tolerance for eigenvalues. */
abstol=0.0;
/* Solve the generalized symmetric eigenvalue problem A*x = lambda*B*x
 * using nag_lapackeig_dsbgvx (f08ubc).
 */
nag_lapackeig_dsbgvx(order,job,Nag_Interval,uplo,n,ka,kb,ab,pdab,bb,
pdbb,q,pdq,vl,vu,il,iu,abstol,&m,w,z,pdz,
index,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from in nag_lapackeig_dsbgvx (f08ubc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Print eigensolution */
printf("Number of eigenvalues found = %8"NAG_IFMT"\n\n",m);
printf(" Eigenvalues\n ");
for(j=0;j<m;++j)
printf(" %7.4f%s",w[j],j%6==5?"\n":" ");
printf("\n");
if(job==Nag_DoBoth){
/* nag_file_print_matrix_real_gen (x04cac): Print Matrix of eigenvectors Z.
 */
printf("\n");
fflush(stdout);
nag_file_print_matrix_real_gen(order,Nag_GeneralMatrix,Nag_NonUnitDiag,n,
m,z,pdz,"Selected eigenvectors",0,
&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n",
fail.message);
exit_status=1;
gotoEND;
}
if(fail.code==NE_CONVERGENCE){
printf("eigenvectors failed to converge\n");
printf("Indices of eigenvectors that did not converge\n ");
for(j=0;j<n;++j)
printf("%8"NAG_IFMT"%s",index[j],j%6==5?"\n":"");
}
}
END:
NAG_FREE(ab);
NAG_FREE(bb);
NAG_FREE(q);
NAG_FREE(w);
NAG_FREE(z);
NAG_FREE(index);
returnexit_status;
}