Sure I'll try,
I was thinking of using ls->B as the A matrix for dgelss, ls->work as
work and ls->Binv as B. The result is then stored in ls->Binv, but in
column-major format.
Right now, the column-major result is transposed in
PetscFVLeastSquaresPseudoInverseSVD_Static, and the row-major result is
copied in the output in PetscFVComputeGradient_LeastSquares. I think
it's because PetscFVLeastSquaresPseudoInverse_Static gives the result in
row-major format.
Would it be alright if I changed
PetscFVLeastSquaresPseudoInverseSVD_Static so that the result would
still be in column-major format ? I could include the result recopy in
the if statement for example.
Moreover, this would be to keep the compatibility with
PetscFVLeastSquaresPseudoInverse_Static, but right now it is manually
disabled (with useSVD = PETSC_TRUE), so I am worrying for nothing ?
Pierre
On 30/09/20 20:38, Jed Brown wrote:
Pierre Seize <[email protected]> writes:
Hi,
In PetscFVLeastSquaresPseudoInverseSVD_Static, there is
Brhs = work;
maxmn = PetscMax(m,n);
for (j=0; j<maxmn; j++) {
for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);
}
where on the calling function, PetscFVComputeGradient_LeastSquares, we
set the arguments m <= numFaces, n <= dim and work <= ls->work. The size
of the work array is computed in PetscFVLeastSquaresSetMaxFaces_LS as:
ls->maxFaces = maxFaces;
m = ls->maxFaces;
n = dim;
nrhs = ls->maxFaces;
minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n),
PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */
It's totally buggy because this formula is for the argument to dgelss, but the
array is being used for a different purpose (to place Brhs).
WORK
WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
LWORK
LWORK is INTEGER
The dimension of the array WORK. LWORK >= 1, and also:
LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )
For good performance, LWORK should generally be larger.
If LWORK = -1, then a workspace query is assumed; the
routine
only calculates the optimal size of the WORK array,
returns
this value as the first entry of the WORK array, and no
error
message related to LWORK is issued by XERBLA.
There should be a separate allocation for Brhs and the work argument should be
passed through to dgelss.
The current code passes
tmpwork = Ainv;
along to dgelss, but we don't know that it's the right size either.
Would you be willing to submit a merge request with your best attempt at fixing
this. I can help review and we'll get it into the 3.14.1 release?
ls->workSize = 5*minwork; /* We can afford to be extra generous */
In my example, the used size (maxmn * maxmn) is 81, and the actual size
(ls->workSize) is 75, and therefore valgrind complains.
Is it because I am missing something, or is it a bug ?
Thanks
Pierre Seize
