f08bpce.c

/* nag_lapackeig_ztpqrt (f08bpc) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<nag.h>
intmain(void){
/* Scalars */
doublernorm;
Integerexit_status=0;
Integerpda,pdb,pdt;
Integeri,j,m,n,nb,nrhs;
/* Arrays */
Complex*a=0,*b=0,*c=0,*t=0;
/* Nag Types */
Nag_OrderTypeorder;
NagErrorfail;
#ifdef NAG_COLUMN_MAJOR
#define A(I, J) a[(J - 1) * pda + I - 1]
#define B(I, J) b[(J - 1) * pdb + I - 1]
#define C(I, J) c[(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]
#define C(I, J) c[(I - 1) * pdb + J - 1]
order=Nag_RowMajor;
#endif
INIT_FAIL(fail);
printf("nag_lapackeig_ztpqrt (f08bpc) Example Program Results\n\n");
fflush(stdout);
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%"NAG_IFMT"%"NAG_IFMT"%"NAG_IFMT"%*[^\n]",&m,&n,&nrhs);
nb=MIN(m,n);
if(!(a=NAG_ALLOC(m*n,Complex))||!(b=NAG_ALLOC(m*nrhs,Complex))||
!(c=NAG_ALLOC(m*nrhs,Complex))||
!(t=NAG_ALLOC(nb*MIN(m,n),Complex))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
#ifdef NAG_COLUMN_MAJOR
pda=m;
pdb=m;
pdt=nb;
#else
pda=n;
pdb=nrhs;
pdt=MIN(m,n);
#endif
/* Read A and 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<=m;++i)
for(j=1;j<=nrhs;++j)
scanf(" ( %lf , %lf )",&B(i,j).re,&B(i,j).im);
scanf("%*[^\n]");
for(i=1;i<=m;++i)
for(j=1;j<=nrhs;++j)
C(i,j)=B(i,j);
/* nag_lapackeig_zgeqrt (f08apc).
 * Compute the QR factorization of first n rows of A by recursive algorithm.
 */
nag_lapackeig_zgeqrt(order,n,n,nb,a,pda,t,pdt,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_zgeqrt (f08apc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_lapackeig_zgemqrt (f08aqc).
 * Compute C = (C1) = (Q^H)*B, storing the result in C
 * (C2)
 * by applying Q^H from left.
 */
nag_lapackeig_zgemqrt(order,Nag_LeftSide,Nag_ConjTrans,n,nrhs,n,nb,a,
pda,t,pdt,c,pdb,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_zgemqrt (f08aqc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
for(i=1;i<=n;++i)
for(j=1;j<=nrhs;++j)
B(i,j)=C(i,j);
/* nag_lapacklin_ztrtrs (f07tsc).
 * Compute least squares solutions for first n rows
 * by back-substitution in R*X = C1.
 */
nag_lapacklin_ztrtrs(order,Nag_Upper,Nag_NoTrans,Nag_NonUnitDiag,n,nrhs,
a,pda,c,pdb,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapacklin_ztrtrs (f07tsc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_file_print_matrix_complex_gen_comp (x04dbc).
 * Print least squares solutions using first n rows.
 */
nag_file_print_matrix_complex_gen_comp(
order,Nag_GeneralMatrix,Nag_NonUnitDiag,n,nrhs,c,pdb,
Nag_BracketForm,"%7.4f","Solution(s) for n rows",Nag_IntegerLabels,0,
Nag_IntegerLabels,0,80,0,0,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_file_print_matrix_complex_gen_comp (x04dbc).\n%s\n",
fail.message);
exit_status=1;
gotoEND;
}
/* nag_lapackeig_ztpqrt (f08bpc).
 * Now add the remaining rows and perform QR update.
 */
nag_lapackeig_ztpqrt(order,m-n,n,0,nb,a,pda,&A(n+1,1),pda,t,
pdt,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_ztpqrt (f08bpc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_lapackeig_ztpmqrt (f08bqc).
 * Apply orthogonal transformations to C.
 */
nag_lapackeig_ztpmqrt(order,Nag_LeftSide,Nag_ConjTrans,m-n,nrhs,n,0,
nb,&A(n+1,1),pda,t,pdt,b,pdb,&B(5,1),pdb,
&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_ztpmqrt (f08bqc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_lapacklin_ztrtrs (f07tsc).
 * Compute least squares solutions for first n rows
 * by back-substitution in R*X = C1.
 */
nag_lapacklin_ztrtrs(order,Nag_Upper,Nag_NoTrans,Nag_NonUnitDiag,n,nrhs,
a,pda,b,pdb,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapacklin_ztrtrs (f07tsc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* nag_file_print_matrix_complex_gen_comp (x04dbc).
 * Print least squares solutions.
 */
printf("\n");
fflush(stdout);
nag_file_print_matrix_complex_gen_comp(
order,Nag_GeneralMatrix,Nag_NonUnitDiag,n,nrhs,b,pdb,
Nag_BracketForm,"%7.4f","Least squares solution(s) for all rows",
Nag_IntegerLabels,0,Nag_IntegerLabels,0,80,0,0,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_file_print_matrix_complex_gen_comp (x04dbc).\n%s\n",
fail.message);
exit_status=1;
gotoEND;
}
printf("\n Square root(s) of the residual sum(s) of squares\n");
for(j=1;j<=nrhs;j++){
/* nag_blast_zge_norm (f16uac).
 * Compute and print estimate of the square root of the residual
 * sum of squares.
 */
nag_blast_zge_norm(order,Nag_FrobeniusNorm,m-n,1,&B(n+1,j),pdb,
&rnorm,&fail);
if(fail.code!=NE_NOERROR){
printf("\nError from nag_blast_zge_norm (f16uac).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
printf(" %11.2e ",rnorm);
}
printf("\n");
END:
NAG_FREE(a);
NAG_FREE(b);
NAG_FREE(c);
NAG_FREE(t);
returnexit_status;
}