f08ysce.c

/* nag_lapackeig_ztgsja (f08ysc) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<nag.h>
#include<stdio.h>
intmain(void){
/* Scalars */
doubleeps,norma,normb,tola,tolb;
Integeri,irank,j,k,l,m,n,ncycle,p,pda,pdb,pdu,pdv;
Integerpdq,printq,printr,printu,printv,vsize;
Integerexit_status=0;
/* Arrays */
Complex*a=0,*b=0,*q=0,*u=0,*v=0;
double*alpha=0,*beta=0;
charnag_enum_arg[40];
/* Nag Types */
NagErrorfail;
Nag_OrderTypeorder;
Nag_ComputeUTypejobu;
Nag_ComputeVTypejobv;
Nag_ComputeQTypejobq;
Nag_MatrixTypegenmat=Nag_GeneralMatrix,upmat=Nag_UpperMatrix;
Nag_DiagTypediag=Nag_NonUnitDiag;
Nag_LabelTypeintlab=Nag_IntegerLabels;
Nag_ComplexFormTypebrac=Nag_BracketForm;
#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_ztgsja (f08ysc) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%"NAG_IFMT"%"NAG_IFMT"%"NAG_IFMT"%*[^\n]",&m,&n,&p);
if(m<0||n<0||p<0){
printf("Invalid m, n or p\n");
exit_status=1;
gotoEND;
}
scanf(" %39s%*[^\n]",nag_enum_arg);
/* nag_enum_name_to_value (x04nac).
 * Converts NAG enum member name to value
 */
jobu=(Nag_ComputeUType)nag_enum_name_to_value(nag_enum_arg);
scanf(" %39s%*[^\n]",nag_enum_arg);
jobv=(Nag_ComputeVType)nag_enum_name_to_value(nag_enum_arg);
scanf(" %39s%*[^\n]",nag_enum_arg);
jobq=(Nag_ComputeQType)nag_enum_name_to_value(nag_enum_arg);
pdu=(jobu!=Nag_NotU?m:1);
pdv=(jobv!=Nag_NotV?p:1);
pdq=(jobq!=Nag_NotQ?n:1);
vsize=(jobv!=Nag_NotV?p*m:1);
#ifdef NAG_COLUMN_MAJOR
pda=m;
pdb=p;
#else
pda=n;
pdb=n;
#endif
/* Read in 0s or 1s to determine whether matrices U, V, Q or R are to be
 * printed.
 */
scanf("%"NAG_IFMT"%"NAG_IFMT"%"NAG_IFMT"%"NAG_IFMT"%*[^\n]",&printu,
&printv,&printq,&printr);
/* Allocate memory */
if(!(a=NAG_ALLOC(m*n,Complex))||!(b=NAG_ALLOC(p*n,Complex))||
!(alpha=NAG_ALLOC(n,double))||!(beta=NAG_ALLOC(n,double))||
!(q=NAG_ALLOC(pdq*pdq,Complex))||
!(u=NAG_ALLOC(pdu*pdu,Complex))||
!(v=NAG_ALLOC(vsize,Complex))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
/* Read the m by n matrix A and p by n matrix B from data file */
for(i=1;i<=m;++i)
for(j=1;j<=n;++j)
scanf(" ( %lf , %lf )",&A(i,j).re,&A(i,j).im);
scanf("%*[^\n]");
for(i=1;i<=p;++i)
for(j=1;j<=n;++j)
scanf(" ( %lf , %lf )",&B(i,j).re,&B(i,j).im);
scanf("%*[^\n]");
/* Compute tola and tolb as */
/* tola = max(m,n)*norm(A)*macheps */
/* tolb = max(p,n)*norm(B)*macheps */
nag_blast_zge_norm(order,Nag_OneNorm,m,n,a,pda,&norma,&fail);
nag_blast_zge_norm(order,Nag_OneNorm,p,n,b,pdb,&normb,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_blast_zge_norm (f16uac).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Compute tola and tolb using nag_machine_precision (x02ajc) */
eps=nag_machine_precision;
tola=MAX(m,n)*norma*eps;
tolb=MAX(p,n)*normb*eps;
/* Preprocess step:
 * compute transformations to reduce (A, B) to upper triangular form
 * (A = U1*S*(Q1^H), B = V1*T*(Q1^H))
 * using nag_lapackeig_zggsvp (f08vsc).
 */
nag_lapackeig_zggsvp(order,jobu,jobv,jobq,m,p,n,a,pda,b,pdb,tola,
tolb,&k,&l,u,pdu,v,pdv,q,pdq,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_zggsvp (f08vsc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Compute the generalized singular value decomposition of preprocessed (A, B)
 * (A = U*D1*(0 R)*(Q^H), B = V*D2*(0 R)*(Q^H))
 * using nag_lapackeig_ztgsja (f08ysc).
 */
nag_lapackeig_ztgsja(order,jobu,jobv,jobq,m,p,n,k,l,a,pda,b,pdb,
tola,tolb,alpha,beta,u,pdu,v,pdv,q,pdq,&ncycle,
&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_ztgsja (f08ysc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Print the generalized singular value pairs alpha, beta */
irank=MIN(k+l,m);
printf("Number of infinite generalized singular values (k): %5"NAG_IFMT"\n",
k);
printf("Number of finite generalized singular values (l): %5"NAG_IFMT"\n",
l);
printf("Effective Numerical rank of (A^H B^HT)^H (k+l): %5"NAG_IFMT"\n",
irank);
printf("\nFinite generalized singular values:\n");
for(j=k;j<irank;++j)
printf("%45s%12.4e\n","",alpha[j]/beta[j]);
printf("\nNumber of cycles of the Kogbetliantz method: %12"NAG_IFMT"\n\n",
ncycle);
if(printu&&jobu!=Nag_NotU){
fflush(stdout);
nag_file_print_matrix_complex_gen_comp(
order,genmat,diag,m,m,u,pdu,brac,"%13.4e","Unitary matrix U",
intlab,NULL,intlab,NULL,80,0,NULL,&fail);
if(fail.code!=NE_NOERROR)
gotoPRINTERR;
}
if(printv&&jobv!=Nag_NotV){
printf("\n");
fflush(stdout);
nag_file_print_matrix_complex_gen_comp(
order,genmat,diag,p,p,v,pdv,brac,"%13.4e","Unitary matrix V",
intlab,NULL,intlab,NULL,80,0,NULL,&fail);
if(fail.code!=NE_NOERROR)
gotoPRINTERR;
}
if(printq&&jobq!=Nag_NotQ){
printf("\n");
fflush(stdout);
nag_file_print_matrix_complex_gen_comp(
order,genmat,diag,n,n,q,pdq,brac,"%13.4e","Unitary matrix Q",
intlab,NULL,intlab,NULL,80,0,NULL,&fail);
if(fail.code!=NE_NOERROR)
gotoPRINTERR;
}
if(printr){
printf("\n");
fflush(stdout);
nag_file_print_matrix_complex_gen_comp(
order,upmat,diag,irank,irank,&A(1,n-irank+1),pda,brac,
"%13.4e","Nonsingular upper triangular matrix R",intlab,NULL,intlab,
NULL,80,0,NULL,&fail);
}
PRINTERR:
if(fail.code!=NE_NOERROR){
printf("Error from nag_file_print_matrix_real_gen_comp (x04cbc).\n%s\n",
fail.message);
exit_status=1;
}
END:
NAG_FREE(a);
NAG_FREE(b);
NAG_FREE(alpha);
NAG_FREE(beta);
NAG_FREE(q);
NAG_FREE(u);
NAG_FREE(v);
returnexit_status;
}