f08bace.c

/* nag_lapackeig_dgelsy (f08bac) Example Program.
 *
 * Copyright 2025 Numerical Algorithms Group.
 *
 * Mark 31.1, 2025.
 */
#include<nag.h>
#include<stdio.h>
intmain(void){
/* Scalars */
doublercond;
Integerexit_status=0,i,j,m,n,nrhs,rank,pda,pdb;
/* Arrays */
double*a=0,*b=0;
Integer*jpvt=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]
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_dgelsy (f08bac) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%"NAG_IFMT"%"NAG_IFMT"%"NAG_IFMT"%*[^\n]",&m,&n,&nrhs);
/* Allocate memory */
if(!(a=NAG_ALLOC(m*n,double))||!(b=NAG_ALLOC(m*nrhs,double))||
!(jpvt=NAG_ALLOC(n,Integer))){
printf("Allocation failure\n");
exit_status=-1;
gotoEND;
}
#ifdef NAG_COLUMN_MAJOR
pda=m;
pdb=MAX(m,n);
#else
pda=n;
pdb=nrhs;
#endif
/* 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]");
/* nag_blast_iload (f16dbc).
 * Initialize jpvt to be zero so that all columns are free.
 */
nag_blast_iload(n,0,jpvt,1,&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_blast_iload (f16dbc).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Choose rcond to reflect the relative accuracy of the input data */
rcond=0.01;
/* nag_lapackeig_dgelsy (f08bac).
 * Solve the least squares problem min( norm2(b - Ax) ) for the x
 * of minimum norm.
 */
nag_lapackeig_dgelsy(order,m,n,nrhs,a,pda,b,pdb,jpvt,rcond,&rank,
&fail);
if(fail.code!=NE_NOERROR){
printf("Error from nag_lapackeig_dgelsy (f08bac).\n%s\n",fail.message);
exit_status=1;
gotoEND;
}
/* Print solution */
printf("Least squares solution\n");
for(i=1;i<=n;++i){
for(j=1;j<=nrhs;++j)
printf("%11.4f%s",B(i,j),j%7==0?"\n":" ");
printf("\n");
}
/* Print the effective rank of A */
printf("\nTolerance used to estimate the rank of A\n");
printf("%11.2e\n",rcond);
printf("Estimated rank of A\n%6"NAG_IFMT"\n",rank);
END:
NAG_FREE(a);
NAG_FREE(b);
NAG_FREE(jpvt);
returnexit_status;
}
#undef A
#undef B