[petsc-users] question about the PetscFVLeastSquaresPseudoInverseSVD
Rongliang Chen
rongliang.chan at gmail.com
Mon Mar 21 22:10:21 CDT 2016
Hello All,
Can anyone help me to understand the output of
PetscFVLeastSquaresPseudoInverseSVD?
I tested six matrices and the results confused me. The results are
followed and the source codes are attached.
For the cases m>n, the results make sense for me, and for the m<=n
cases, I do not understand the results. For example, for the m=n=3
cases, I think I should get an identity matrix from this function and I
get two difference outputs for the same input matrix.
Best regards,
Rongliang
-----------------------------------
Initialize the matrix A (is a 3 x 4 matrix):
1.000000 0.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
Initialize the matrix B:
1.000000 0.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
0.000000 0.000000 0.000000 1.000000
The output of the SVD based least square:
1.000000 0.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
-----------------------------------
Initialize the matrix A (is a 3 x 4 matrix):
1.000000 0.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
Initialize the matrix B:
1.000000 0.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
0.000000 0.000000 0.000000 1.000000
The output of the SVD based least square:
1.000000 0.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
-----------------------------------
Initialize the matrix A (is a 3 x 3 matrix):
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
Initialize the matrix B:
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
The output of the SVD based least square:
0.500000 0.000000 -0.500000
0.500000 0.000000 -0.500000
0.000000 0.000000 1.000000
-----------------------------------
Initialize the matrix A (is a 3 x 3 matrix):
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
Initialize the matrix B:
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
The output of the SVD based least square:
0.500000 -0.707107 0.000000
0.500000 -0.707107 0.000000
-0.000000 1.414214 0.000000
-----------------------------------
Initialize the matrix A (is a 3 x 2 matrix):
1.000000 0.000000
0.000000 1.000000
0.000000 0.000000
Initialize the matrix B:
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
The output of the SVD based least square:
1.000000 -0.838728
0.000000 0.000000
-0.000000 1.305169
-----------------------------------
Initialize the matrix A (is a 3 x 2 matrix):
1.000000 0.000000
0.000000 1.000000
0.000000 0.000000
Initialize the matrix B:
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
The output of the SVD based least square:
1.000000 -0.293611
0.000000 0.000000
-0.000000 1.000000
-------------- next part --------------
#include <petscts.h>
#include <petscfv.h>
#include <petscdmplex.h>
#include <petscsf.h>
#include <petscblaslapack.h>
#include <petsctime.h>
static char help[] = "Test the SVD function.\n";
#undef __FUNCT__
#define __FUNCT__ "PetscFVLeastSquaresPseudoInverseSVD"
/* Overwrites A. Can handle degenerate problems and m<n. */
PetscErrorCode PetscFVLeastSquaresPseudoInverseSVD(PetscInt m,PetscInt mstride,PetscInt n,PetscScalar *A,PetscScalar *Ainv,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
{
PetscBool debug = PETSC_FALSE;
PetscScalar *Brhs, *Aback;
PetscScalar *tmpwork;
PetscReal rcond;
PetscInt i, j, maxmn;
PetscBLASInt M, N, lda, ldb, ldwork;
#if !defined(PETSC_USE_COMPLEX)
PetscBLASInt nrhs, irank, info;
#endif
PetscErrorCode ierr;
PetscMPIInt rank;
PetscFunctionBegin;
ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
if (debug) {
ierr = PetscMalloc1(m*n,&Aback);CHKERRQ(ierr);
ierr = PetscMemcpy(Aback,A,m*n*sizeof(PetscScalar));CHKERRQ(ierr);
}
tmpwork = Ainv;
Brhs = work;
maxmn = PetscMax(m,n);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Initialize the matrix A (is a %d x %d matrix):\n", n, m);
/* initialize to identity */
for (j=0; j<n; j++) {
for (i=0; i<m; i++) {
if (i == j) {
A[i + j*m] = 1.0;
}else{
A[i + j*m] = 0.0;
}
PetscPrintf(PETSC_COMM_WORLD,"%f ", A[i + j*m]);
}
PetscPrintf(PETSC_COMM_WORLD,"\n");
}
PetscPrintf(PETSC_COMM_WORLD,"\n");
ierr = PetscPrintf(PETSC_COMM_WORLD,"Initialize the matrix B:\n");
/* initialize to identity */
if (rank==2) PetscPrintf(PETSC_COMM_WORLD,"m = %d, n = %d, maxmn = %d\n", m, n, maxmn);
for (j=0; j<maxmn; j++) {
for (i=0; i<maxmn; i++) {
if (i == j) {
Brhs[i + j*maxmn] = 1.0;
}else{
Brhs[i + j*maxmn] = 0.0;
}
PetscPrintf(PETSC_COMM_WORLD,"%f ", Brhs[i + j*maxmn]);
}
PetscPrintf(PETSC_COMM_WORLD,"\n");
}
PetscPrintf(PETSC_COMM_WORLD,"\n");
ierr = PetscBLASIntCast(m,&M);CHKERRQ(ierr);
ierr = PetscBLASIntCast(n,&N);CHKERRQ(ierr);
ierr = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr);
ierr = PetscBLASIntCast(maxmn,&ldb);CHKERRQ(ierr);
ierr = PetscBLASIntCast(worksize,&ldwork);CHKERRQ(ierr);
rcond = -1;
ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
#if defined(PETSC_USE_COMPLEX)
if (tmpwork && rcond) rcond = 0.0; /* Get rid of compiler warning */
SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "I don't think this makes sense for complex numbers");
#else
nrhs = M;
LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);
ierr = PetscFPTrapPop();CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"The output of the SVD based least square:\n");
for (i=0; i<n; i++) {
for (j=0; j<nrhs; j++) {
Ainv[i*mstride+j] = Brhs[i + j*maxmn];
PetscPrintf(PETSC_COMM_WORLD,"%f ", Ainv[i*mstride+j]);
}
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");
}
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n");
PetscFunctionReturn(0);
#endif
}
#undef __FUNCT__
#define __FUNCT__ "PseudoInverseGetWorkRequired"
PetscErrorCode PseudoInverseGetWorkRequired(PetscInt maxFaces,PetscInt *work)
{
PetscInt m,n,nrhs,minwork;
PetscFunctionBeginUser;
m = maxFaces;
n = 3;
nrhs = maxFaces;
minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n), PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
*work = 5*minwork; /* We can afford to be extra generous */
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv)
{
PetscErrorCode ierr;
PetscInt maxNumFaces = 4, worksize;
PetscScalar *B,*Binv,*work,*tau;
PetscInt dim = 3, numFaces = 4;
PetscInt usedFaces[6] = {4, 4, 3, 3, 2, 2};
PetscInt i;
ierr = PetscInitialize(&argc, &argv, (char*) 0, help);CHKERRQ(ierr);
ierr = PseudoInverseGetWorkRequired(maxNumFaces,&worksize);CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD, "maxNumFaces = %d, worksize = %d\n", maxNumFaces, worksize);
ierr = PetscMalloc4(maxNumFaces*dim, &B, worksize, &Binv, worksize, &work, maxNumFaces, &tau);CHKERRQ(ierr);
for(i=0; i<6; i++){
PetscPrintf(PETSC_COMM_WORLD, "-----------------------------------\n");
ierr = PetscFVLeastSquaresPseudoInverseSVD(usedFaces[i],numFaces,dim,B,Binv,tau,worksize,work);CHKERRQ(ierr);
}
ierr = PetscFree4(B,Binv,work,tau);CHKERRQ(ierr);
ierr = PetscFinalize();
return(0);
}
More information about the petsc-users
mailing list