<div dir="ltr"><div dir="ltr">On Thu, Oct 1, 2020 at 4:18 AM Pierre Seize <<a href="mailto:Pierre.Seize@onera.fr">Pierre.Seize@onera.fr</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Sure I'll try,<br>
<br>
I was thinking of using ls->B as the A matrix for dgelss, ls->work as <br>
work and ls->Binv as B. The result is then stored in ls->Binv, but in <br>
column-major format.<br>
<br>
Right now, the column-major result is transposed in <br>
PetscFVLeastSquaresPseudoInverseSVD_Static, and the row-major result is <br>
copied in the output in PetscFVComputeGradient_LeastSquares. I think <br>
it's because PetscFVLeastSquaresPseudoInverse_Static gives the result in <br>
row-major format.<br>
<br>
Would it be alright if I changed <br>
PetscFVLeastSquaresPseudoInverseSVD_Static so that the result would <br>
still be in column-major format ? I could include the result recopy in <br>
the if statement for example.<br>
<br>
Moreover, this would be to keep the compatibility with <br>
PetscFVLeastSquaresPseudoInverse_Static, but right now it is manually <br>
disabled (with useSVD = PETSC_TRUE), so I am worrying for nothing ?<br></blockquote><div><br></div><div>As long as the layout is documented in the function manpage, I think that change is fine. Nothing else uses</div><div>the code except the reconstruction right now.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Pierre<br>
<br>
<br>
On 30/09/20 20:38, Jed Brown wrote:<br>
> Pierre Seize <<a href="mailto:Pierre.Seize@onera.fr" target="_blank">Pierre.Seize@onera.fr</a>> writes:<br>
><br>
>> Hi,<br>
>><br>
>> In PetscFVLeastSquaresPseudoInverseSVD_Static, there is<br>
>> Brhs = work;<br>
>> maxmn = PetscMax(m,n);<br>
>> for (j=0; j<maxmn; j++) {<br>
>> for (i=0; i<maxmn; i++) Brhs[i + j*maxmn] = 1.0*(i == j);<br>
>> }<br>
>> where on the calling function, PetscFVComputeGradient_LeastSquares, we<br>
>> set the arguments m <= numFaces, n <= dim and work <= ls->work. The size<br>
>> of the work array is computed in PetscFVLeastSquaresSetMaxFaces_LS as:<br>
>> ls->maxFaces = maxFaces;<br>
>> m = ls->maxFaces;<br>
>> n = dim;<br>
>> nrhs = ls->maxFaces;<br>
>> minwork = 3*PetscMin(m,n) + PetscMax(2*PetscMin(m,n),<br>
>> PetscMax(PetscMax(m,n), nrhs)); /* required by LAPACK */<br>
> 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).<br>
><br>
> WORK<br>
><br>
> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))<br>
> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.<br>
><br>
> LWORK<br>
><br>
> LWORK is INTEGER<br>
> The dimension of the array WORK. LWORK >= 1, and also:<br>
> LWORK >= 3*min(M,N) + max( 2*min(M,N), max(M,N), NRHS )<br>
> For good performance, LWORK should generally be larger.<br>
><br>
> If LWORK = -1, then a workspace query is assumed; the routine<br>
> only calculates the optimal size of the WORK array, returns<br>
> this value as the first entry of the WORK array, and no error<br>
> message related to LWORK is issued by XERBLA.<br>
><br>
> There should be a separate allocation for Brhs and the work argument should be passed through to dgelss.<br>
><br>
> The current code passes<br>
><br>
> tmpwork = Ainv;<br>
><br>
> along to dgelss, but we don't know that it's the right size either.<br>
><br>
><br>
> 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?<br>
><br>
>> ls->workSize = 5*minwork; /* We can afford to be extra generous */<br>
>><br>
>> In my example, the used size (maxmn * maxmn) is 81, and the actual size<br>
>> (ls->workSize) is 75, and therefore valgrind complains.<br>
>> Is it because I am missing something, or is it a bug ?<br>
>><br>
>> Thanks<br>
>><br>
>> Pierre Seize<br>
<br>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>