f08bece.c

/* nag_lapackeig_dgeqpf (f08bec) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<nag.h>
#include<stdio.h>
intmain(void){
/* Scalars */
doubletol;
Integeri,j,jpvt_len,k,m,n,nrhs;
Integerpda,pdb,pdx,tau_len;
Integerexit_status=0;
NagErrorfail;
Nag_OrderTypeorder;
/* Arrays */
double*a=0,*b=0,*tau=0,*x=0;
Integer*jpvt=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]
#define X(I, J) x[(J - 1) * pdx + 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 X(I, J) x[(I - 1) * pdx + J - 1]
order=Nag_RowMajor;
#endif
INIT_FAIL(fail);
printf("nag_lapackeig_dgeqpf (f08bec) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n] ");
scanf("%"NAG_IFMT"%"NAG_IFMT"%"NAG_IFMT"%*[^\n] ",&m,&n,&nrhs);
#ifdef NAG_COLUMN_MAJOR
pda=m;
pdb=m;
pdx=m;
#else
pda=n;
pdb=nrhs;
pdx=nrhs;
#endif
tau_len=MIN(m,n);
jpvt_len=n;
/* Allocate memory */
if(!(a=NAG_ALLOC(m*n,double))||!(b=NAG_ALLOC(m*nrhs,double))||
!(tau=NAG_ALLOC(tau_len,double))||
!(x=NAG_ALLOC(m*nrhs,double))||
!(jpvt=NAG_ALLOC(jpvt_len,Integer))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
/* Read A and B from data file */
for(i=1;i<=m;++i){
for(j=1;j<=n;++j)
scanf("%lf",&A(i,j));
}
scanf("%*[^\n] ");
for(i=1;i<=m;++i){
for(j=1;j<=nrhs;++j)
scanf("%lf",&B(i,j));
}
scanf("%*[^\n] ");
/* Initialize JPVT to be zero so that all columns are free */
/* nag_blast_iload (f16dbc).
 * Broadcast scalar into integer vector
 */
nag_blast_iload(n,0,jpvt,1,&fail);
/* Compute the QR factorization of A */
/* nag_lapackeig_dgeqpf (f08bec).
 * QR factorization of real general rectangular matrix with
 * column pivoting
 */
nag_lapackeig_dgeqpf(order,m,n,a,pda,jpvt,tau,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dgeqpf (f08bec).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Choose TOL to reflect the relative accuracy of the input data */
tol=0.01;
/* Determine which columns of R to use */
for(k=1;k<=n;++k){
if(ABS(A(k,k))<=tol*ABS(A(1,1)))
break;
}
--k;
/* Compute C = (Q^T)*B, storing the result in B */
/* 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,m,nrhs,n,a,pda,tau,
b,pdb,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dormqr (f08agc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Compute least squares solution by back-substitution in R*B = C */
/* nag_lapacklin_dtrtrs (f07tec).
 * Solution of real triangular system of linear equations,
 * multiple right-hand sides
 */
nag_lapacklin_dtrtrs(order,Nag_Upper,Nag_NoTrans,Nag_NonUnitDiag,k,nrhs,
a,pda,b,pdb,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapacklin_dtrtrs (f07tec).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
for(i=k+1;i<=n;++i){
for(j=1;j<=nrhs;++j)
B(i,j)=0.0;
}
/* Unscramble the least squares solution stored in B */
for(i=1;i<=n;++i){
for(j=1;j<=nrhs;++j)
X(jpvt[i-1],j)=B(i,j);
}
/* Print least squares solution */
/* 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,
nrhs,x,pdx,"Least squares solution",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(b);
NAG_FREE(tau);
NAG_FREE(x);
NAG_FREE(jpvt);
returnexit_status;
}