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
