[petsc-users] question about the PetscFVLeastSquaresPseudoInverseSVD

[Jed Brown]

Rongliang Chen <rongliang.chan at gmail.com> writes:

> -----------------------------------
> 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

So the above looks like a problem.  Let's see how that matrix is
constructed.

>   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]);

It's a packed 3x3 matrix (lda=3).

>   ierr  = PetscBLASIntCast(mstride,&lda);CHKERRQ(ierr);

And yet mstride=4, so your matrix is not packed correctly.  If you're
not familiar with BLAS-style packing with lda, please read the
documentation.

>   LAPACKgelss_(&M,&N,&nrhs,A,&lda,Brhs,&ldb, (PetscReal *) tau,&rcond,&irank,tmpwork,&ldwork,&info);

If you use A[i + j*mstride] when filling in the entries above, the
program outputs the following.

maxNumFaces = 4, worksize = 75
-----------------------------------
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 0.000000 1.000000 
0.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.000000 0.000000 
0.000000 1.000000 0.000000 
0.000000 0.000000 1.000000 

-----------------------------------
Initialize the matrix A (is a 3 x 3 matrix):
1.000000 0.000000 0.000000 
0.000000 0.000000 1.000000 
0.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.000000 0.000000 
0.000000 1.000000 0.000000 
0.000000 0.000000 1.000000 

-----------------------------------
Initialize the matrix A (is a 3 x 2 matrix):
1.000000 0.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:
1.000000 0.000000 
0.000000 1.000000 
0.000000 0.000000 

-----------------------------------
Initialize the matrix A (is a 3 x 2 matrix):
1.000000 0.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:
1.000000 0.000000 
0.000000 1.000000 
0.000000 0.000000 

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 818 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160329/b1fadef7/attachment.pgp>