f08nhce.c

/* nag_lapackeig_dgebal (f08nhc) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<nag.h>
#include<stdio.h>
intmain(void){
/* Scalars */
Integerfirstnz,i,ihi,ilo,j,m,n,pda,pdh,pdvr;
Integerscale_len,tau_len,wi_len,wr_len;
Integerexit_status=0;
NagErrorfail;
Nag_OrderTypeorder;
/* Arrays */
double*a=0,*h=0,*scale=0,*tau=0,*vl=0,*vr=0;
double*wi=0,*wr=0;
Nag_Boolean*select=0;
#ifdef NAG_COLUMN_MAJOR
#define A(I, J) a[(J - 1) * pda + I - 1]
#define H(I, J) h[(J - 1) * pdh + 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 H(I, J) h[(I - 1) * pdh + J - 1]
#define VR(I, J) vr[(I - 1) * pdvr + J - 1]
order=Nag_RowMajor;
#endif
INIT_FAIL(fail);
printf("nag_lapackeig_dgebal (f08nhc) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%"NAG_IFMT"%*[^\n] ",&n);
pda=n;
pdh=n;
pdvr=n;
scale_len=n;
tau_len=n;
wi_len=n;
wr_len=n;
/* Allocate memory */
if(!(a=NAG_ALLOC(n*n,double))||!(h=NAG_ALLOC(n*n,double))||
!(scale=NAG_ALLOC(scale_len,double))||
!(tau=NAG_ALLOC(tau_len,double))||!(vl=NAG_ALLOC(1*1,double))||
!(vr=NAG_ALLOC(n*n,double))||!(wi=NAG_ALLOC(wi_len,double))||
!(wr=NAG_ALLOC(wr_len,double))||
!(select=NAG_ALLOC(1,Nag_Boolean))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
/* Read A from data file */
for(i=1;i<=n;++i){
for(j=1;j<=n;++j)
scanf("%lf",&A(i,j));
}
scanf("%*[^\n] ");
/* Balance A */
/* nag_lapackeig_dgebal (f08nhc).
 * Balance real general matrix
 */
nag_lapackeig_dgebal(order,Nag_DoBoth,n,a,pda,&ilo,&ihi,scale,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dgebal (f08nhc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Reduce A to upper Hessenberg form H = (Q^T)*A*Q */
/* nag_lapackeig_dgehrd (f08nec).
 * Orthogonal reduction of real general matrix to upper
 * Hessenberg form
 */
nag_lapackeig_dgehrd(order,n,ilo,ihi,a,pda,tau,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dgehrd (f08nec).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Copy A to H and VR */
for(i=1;i<=n;++i){
for(j=1;j<=n;++j){
H(i,j)=A(i,j);
VR(i,j)=A(i,j);
}
}
/* Form Q explicitly, storing the result in VR */
/* nag_lapackeig_dorghr (f08nfc).
 * Generate orthogonal transformation matrix from reduction
 * to Hessenberg form determined by nag_lapackeig_dgehrd (f08nec)
 */
nag_lapackeig_dorghr(order,n,1,n,vr,pdvr,tau,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dorghr (f08nfc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Calculate the eigenvalues and Schur factorization of A */
/* nag_lapackeig_dhseqr (f08pec).
 * Eigenvalues and Schur factorization of real upper
 * Hessenberg matrix reduced from real general matrix
 */
nag_lapackeig_dhseqr(order,Nag_Schur,Nag_UpdateZ,n,ilo,ihi,h,pdh,wr,
wi,vr,pdvr,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dhseqr (f08pec).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
printf(" Eigenvalues\n");
for(i=1;i<=n;++i)
printf("(%8.4f,%8.4f)\n",wr[i-1],wi[i-1]);
/* Calculate the eigenvectors of A, storing the result in VR */
/* nag_lapackeig_dtrevc (f08qkc).
 * Left and right eigenvectors of real upper
 * quasi-triangular matrix
 */
nag_lapackeig_dtrevc(order,Nag_RightSide,Nag_BackTransform,select,n,h,
pdh,vl,1,vr,pdvr,n,&m,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dtrevc (f08qkc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_lapackeig_dgebak (f08njc).
 * Transform eigenvectors of real balanced matrix to those
 * of original matrix supplied to nag_lapackeig_dgebal (f08nhc)
 */
nag_lapackeig_dgebak(order,Nag_DoBoth,Nag_RightSide,n,ilo,ihi,scale,m,
vr,pdvr,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dgebak (f08njc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Normalize the left eigenvectors */
for(j=1;j<=n;j++){
firstnz=n;
for(i=n;i>=1;i--){
if(VR(i,j)!=0){
firstnz=i;
}
}
for(i=n;i>=1;i--){
VR(i,j)=VR(i,j)/VR(firstnz,j);
}
}
/* Print eigenvectors */
printf("\n");
/* nag_file_print_matrix_real_gen (x04cac).
 * Print real general matrix (easy-to-use)
 */
fflush(stdout);
nag_file_print_matrix_real_gen(order,Nag_GeneralMatrix,Nag_NonUnitDiag,n,
m,vr,pdvr,"Contents of array VR",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;
}
END:
NAG_FREE(a);
NAG_FREE(h);
NAG_FREE(scale);
NAG_FREE(tau);
NAG_FREE(vl);
NAG_FREE(vr);
NAG_FREE(wi);
NAG_FREE(wr);
NAG_FREE(select);
returnexit_status;
}