f08fnce.c

/* nag_lapackeig_zheev (f08fnc) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<math.h>
#include<nag.h>
#include<stdio.h>
intmain(void){
/* Scalars */
doubleeerrbd,eps;
Integerexit_status=0,i,j,n,pda;
/* Arrays */
Complex*a=0;
double*rcondz=0,*w=0,*zerrbd=0;
/* Nag Types */
Nag_OrderTypeorder;
NagErrorfail;
#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_zheev (f08fnc) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%"NAG_IFMT"%*[^\n]",&n);
/* Allocate memory */
if(!(a=NAG_ALLOC(n*n,Complex))||!(rcondz=NAG_ALLOC(n,double))||
!(w=NAG_ALLOC(n,double))||!(zerrbd=NAG_ALLOC(n,double))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
#ifdef NAG_COLUMN_MAJOR
pda=n;
#else
pda=n;
#endif
/* Read the upper triangular part of the matrix A from data file */
for(i=1;i<=n;++i)
for(j=i;j<=n;++j)
scanf(" ( %lf , %lf )",&A(i,j).re,&A(i,j).im);
scanf("%*[^\n]");
/* nag_lapackeig_zheev (f08fnc).
 * Solve the Hermitian eigenvalue problem.
 */
nag_lapackeig_zheev(order,Nag_DoBoth,Nag_Upper,n,a,pda,w,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_zheev (f08fnc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_complex_divide (a02cdc).
 * Normalize the eigenvectors.
 */
for(j=1;j<=n;j++)
for(i=n;i>=1;i--)
A(i,j)=nag_complex_divide(A(i,j),A(1,j));
/* Print solution */
printf("Eigenvalues\n");
for(j=0;j<n;++j)
printf("%8.4f%s",w[j],(j+1)%8==0?"\n":" ");
printf("\n\n");
/* nag_file_print_matrix_complex_gen (x04dac).
 * Print eigenvectors.
 */
fflush(stdout);
nag_file_print_matrix_complex_gen(order,Nag_GeneralMatrix,Nag_NonUnitDiag,
n,n,a,pda,"Eigenvectors",0,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_file_print_matrix_complex_gen (x04dac).\n%s\n",
fail.message);
exit_status=1;
gotoEND;
}
/* Get the machine precision, eps, using nag_machine_precision (X02AJC)
 * and compute the approximate error bound for the computed eigenvalues.
 * Note that for the 2-norm, ||A|| = max {|w[i]|, i=0..n-1}, and since
 * the eigenvalues are in ascending order ||A|| = max( |w[0]|, |w[n-1]|).
 */
eps=nag_machine_precision;
eerrbd=eps*MAX(fabs(w[0]),fabs(w[n-1]));
/* nag_lapackeig_ddisna (f08flc).
 * Estimate reciprocal condition numbers for the eigenvectors.
 */
nag_lapackeig_ddisna(Nag_EigVecs,n,n,w,rcondz,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_ddisna (f08flc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Compute the error estimates for the eigenvectors */
for(i=0;i<n;++i)
zerrbd[i]=eerrbd/rcondz[i];
/* Print the approximate error bounds for the eigenvalues and vectors */
printf("\nError estimate for the eigenvalues\n");
printf("%11.1e\n\n",eerrbd);
printf("Error estimates for the eigenvectors\n");
for(i=0;i<n;++i)
printf("%11.1e%s",zerrbd[i],(i+1)%6==0||i==n-1?"\n":" ");
END:
NAG_FREE(a);
NAG_FREE(rcondz);
NAG_FREE(w);
NAG_FREE(zerrbd);
returnexit_status;
}
#undef A