f08xece.c

/* nag_lapackeig_dhgeqz (f08xec) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<nag.h>
#include<stdio.h>
intmain(void){
/* Scalars */
Integeri,ihi,ilo,inca,irows,j,jlen,k,n,pda,pdb;
Integeralpha_len,beta_len,scale_len,tau_len;
Integerexit_status=0;
doubler;
NagErrorfail;
Nag_OrderTypeorder;
/* Arrays */
double*a=0,*alphai=0,*alphar=0,*b=0,*beta=0;
double*lscale=0,*q=0,*rscale=0,*tau=0,*z=0;
#ifdef NAG_COLUMN_MAJOR
#define A(I, J) a[(J - 1) * pda + I - 1]
#define B(I, J) b[(J - 1) * pdb + 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]
order=Nag_RowMajor;
#endif
INIT_FAIL(fail);
printf("nag_lapackeig_dhgeqz (f08xec) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%"NAG_IFMT"%*[^\n] ",&n);
pda=n;
pdb=n;
alpha_len=n;
beta_len=n;
scale_len=n;
tau_len=n;
#ifdef NAG_COLUMN_MAJOR
inca=1;
#else
inca=pda;
#endif
/* Allocate memory */
if(!(a=NAG_ALLOC(n*n,double))||
!(alphai=NAG_ALLOC(alpha_len,double))||
!(alphar=NAG_ALLOC(alpha_len,double))||
!(b=NAG_ALLOC(n*n,double))||
!(beta=NAG_ALLOC(beta_len,double))||
!(lscale=NAG_ALLOC(scale_len,double))||
!(q=NAG_ALLOC(1*1,double))||
!(rscale=NAG_ALLOC(scale_len,double))||
!(tau=NAG_ALLOC(tau_len,double))||!(z=NAG_ALLOC(1*1,double))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
/* READ matrix A from data file */
for(i=1;i<=n;++i){
for(j=1;j<=n;++j)
scanf("%lf",&A(i,j));
}
scanf("%*[^\n] ");
/* READ matrix B from data file */
for(i=1;i<=n;++i){
for(j=1;j<=n;++j)
scanf("%lf",&B(i,j));
}
scanf("%*[^\n] ");
/* Balance matrix pair (A,B) */
/* nag_lapackeig_dggbal (f08whc).
 * Balance a pair of real general matrices
 */
nag_lapackeig_dggbal(order,Nag_DoBoth,n,a,pda,b,pdb,&ilo,&ihi,lscale,
rscale,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dggbal (f08whc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Matrix A after balancing */
/* 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,
n,a,pda,"Matrix A after balancing",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;
}
printf("\n");
/* Matrix B after balancing */
/* nag_file_print_matrix_real_gen (x04cac), see above. */
fflush(stdout);
nag_file_print_matrix_real_gen(order,Nag_GeneralMatrix,Nag_NonUnitDiag,n,
n,b,pdb,"Matrix B after balancing",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;
}
printf("\n");
/* Reduce B to triangular form using QR */
irows=ihi+1-ilo;
/* nag_lapackeig_dgeqrf (f08aec).
 * QR factorization of real general rectangular matrix
 */
nag_lapackeig_dgeqrf(order,irows,irows,&B(ilo,ilo),pdb,tau,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dgeqrf (f08aec).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Apply the orthogonal transformation to matrix A */
/* nag_lapackeig_dormqr (f08agc).
 * Apply orthogonal transformation determined by nag_lapackeig_dgeqrf
 * (f08aec) or nag_lapackeig_dgeqpf (f08bec)
 */
nag_lapackeig_dormqr(order,Nag_LeftSide,Nag_Trans,irows,irows,irows,
&B(ilo,ilo),pdb,tau,&A(ilo,ilo),pda,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dormqr (f08agc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Compute the generalized Hessenberg form of (A,B) */
/* nag_lapackeig_dgghd3 (f08wfc).
 * Orthogonal reduction of a pair of real general matrices
 * to generalized upper Hessenberg form
 */
nag_lapackeig_dgghd3(order,Nag_NotQ,Nag_NotZ,irows,1,irows,&A(ilo,ilo),
pda,&B(ilo,ilo),pdb,q,1,z,1,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dgghd3 (f08wfc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* A and B can have corresonding columns negated. To standardize, select sign
 * of corresponding columns so that largest element of A column is positive.
 */
for(j=1;j<=n;++j){
jlen=j+1;
if(jlen>n){
jlen=n;
}
nag_blast_damax_val(jlen,&A(1,j),inca,&k,&r,&fail);
if(A(k+1,j)<0.0){
for(i=1;i<=jlen;++i){
A(i,j)=-A(i,j);
}
for(i=1;i<=j;++i){
B(i,j)=-B(i,j);
}
}
}
/* Matrix A in generalized Hessenberg form */
/* nag_file_print_matrix_real_gen (x04cac), see above. */
fflush(stdout);
nag_file_print_matrix_real_gen(order,Nag_GeneralMatrix,Nag_NonUnitDiag,n,
n,a,pda,"Matrix A in Hessenberg form",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;
}
printf("\n");
/* Matrix B in generalized Hessenberg form */
/* nag_file_print_matrix_real_gen (x04cac), see above. */
fflush(stdout);
nag_file_print_matrix_real_gen(order,Nag_GeneralMatrix,Nag_NonUnitDiag,n,
n,b,pdb,"Matrix B is triangular",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;
}
/* Compute the generalized Schur form */
/* nag_lapackeig_dhgeqz (f08xec).
 * Eigenvalues and generalized Schur factorization of real
 * generalized upper Hessenberg form reduced from a pair of
 * real general matrices
 */
nag_lapackeig_dhgeqz(order,Nag_EigVals,Nag_NotQ,Nag_NotZ,n,ilo,ihi,a,
pda,b,pdb,alphar,alphai,beta,q,1,z,1,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dhgeqz (f08xec).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Print the generalized eigenvalues */
printf("\n Generalized eigenvalues\n");
for(i=1;i<=n;++i){
if(beta[i-1]!=0.0){
printf(" %4"NAG_IFMT" (%7.3f,%7.3f)\n",i,
alphar[i-1]/beta[i-1],alphai[i-1]/beta[i-1]);
}else
printf(" %4"NAG_IFMT"Eigenvalue is infinite\n",i);
}
END:
NAG_FREE(a);
NAG_FREE(alphai);
NAG_FREE(alphar);
NAG_FREE(b);
NAG_FREE(beta);
NAG_FREE(lscale);
NAG_FREE(q);
NAG_FREE(rscale);
NAG_FREE(tau);
NAG_FREE(z);
returnexit_status;
}