[petsc-users] Memory violation in PetscFVLeastSquaresPseudoInverseSVD_Static

[Pierre Seize]

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 <Pierre.Seize at onera.fr> 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