[petsc-users] Problem setting Jacobian matrix in complementarity problem
Christopher Athaide
cnathaide at yahoo.com
Wed Apr 13 13:38:15 CDT 2016
Barry,
Thanks. I was able to get the MatCopy to work inside the Jacobian and
produce the optimal results.
I initially created a copy of the Jacobian matrix element-by-element and
stored a copy in the AppCtx data structure. In this way, the diagonal
elements were saved even if they were zero.
There was drastic improvement in the running time. In the original run, the
total run time was 73 seconds, and 57 seconds of these was spent in the
FormJacobian function. Using a MatCopy, the total run time was about 18
seconds.
Thanks again for your assistance.
Chris
-----Original Message-----
From: Barry Smith [mailto:bsmith at mcs.anl.gov]
Sent: Tuesday, April 12, 2016 1:08 PM
To: Christopher Athaide <cnathaide at yahoo.com>
Cc: petsc-users at mcs.anl.gov
Subject: Re: [petsc-users] Problem setting Jacobian matrix in
complementarity problem
Yes call the MatCopy() inside the FormJacobian
> On Apr 12, 2016, at 9:32 AM, Christopher Athaide <cnathaide at yahoo.com>
wrote:
>
> Barry,
>
> Thank you for your suggestions. I implemented the changes that you
> suggested.
>
> int main()
> {
> PetscScalar v;
> ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 5 * N, NULL, &J);
> CHKERRQ(ierr);
> for (auto i = 0; i != N; ++i)
> {
> for (auto j = 0; j != N; ++j)
> {
> ierr = MatGetValue(data.A, i, j, &v); CHKERRQ(ierr);
> ierr = MatSetValue(J, i, j, v, INSERT_VALUES);
CHKERRQ(ierr);
> }
> }
>
> /* Assemble the matrix */
> ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> compare_matrices(data.A, J);
>
> }
>
> PetscErrorCode FormJacobian(Tao, Vec, Mat H, Mat, void* ptr) {
>
> PetscFunctionBegin;
>
> // AppCtx *data = (AppCtx *)ptr;
> // Is there any way I can do a copy data->A to H using a matrix or a
> batch operation?
>
> PetscFunctionReturn(0);
> }
>
>
> The good news is that the model runs without generating any errors.
> The bad news is that the program terminates without the optimality
> conditions being satisfied.
>
> In the original program, the model takes about 70 seconds to run and
> the optimality conditions are satisfied. In this version, it takes
> about 5 seconds.
>
> Is there any way I can copy the matrix data->A to H in the
> FormJacobian using a matrix function instead of iterating through all
> the rows and columns and then assembling the matrix at every iteration.
>
> Thank you for your help.
>
> Chris
>
> -----Original Message-----
> From: Barry Smith [mailto:bsmith at mcs.anl.gov]
> Sent: Monday, April 11, 2016 4:05 PM
> To: Christopher Athaide <cnathaide at yahoo.com>
> Cc: petsc-users at mcs.anl.gov
> Subject: Re: [petsc-users] Problem setting Jacobian matrix in
> complementarity problem
>
>
> Does your matrix have entries all along the diagonal? If not you
> need to preallocate for even the 0 locations and set the value of 0 in
> there specifically.
>
> Barry
>
>> On Apr 11, 2016, at 3:59 PM, Christopher Athaide
>> <cnathaide at yahoo.com>
> wrote:
>>
>> Hi Barry,
>>
>> Thank you for taking a look at the issue I am having. I followed your
>> suggestion and made the following changes:
>>
>> main()
>> {
>> ierr = MatDuplicate(data.A, MAT_SHARE_NONZERO_PATTERN, &J);
>> CHKERRQ(ierr);
>> ierr = MatCopy(data.A, J, SAME_NONZERO_PATTERN); CHKERRQ(ierr);
>>
>> }
>>
>> PetscErrorCode FormJacobian(Tao, Vec, Mat, Mat, void*) {
>>
>> PetscFunctionBegin;
>>
>> PetscFunctionReturn(0);
>> }
>>
>> But I get a bunch of errors - I am not sure why.
>>
>> Problem size : 259 rows
>> [0]PETSC ERROR : -------------------- - Error
>> Message--------------------------------------------------------------
>> [0]PETSC ERROR : Argument out of range [0]PETSC ERROR : New nonzero
>> at(0, 0) caused a malloc Use MatSetOption(A,
>> MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check
>> [0]PETSC ERROR : See http
>> ://www.mcs.anl.gov/petsc/documentation/faq.html
>> for trouble shooting.
>> [0]PETSC ERROR : Petsc Release Version 3.6.3, Dec, 03, 2015 [0]PETSC
>> ERROR : Unknown Name on a arch - mswin - c - debug named LAGRANGE by
>> Christopher Mon Apr 11 21 : 50 : 06 2016 [0]PETSC ERROR : Configure
>> options --with - cc = "win32fe cl" --with - fc =
>> 0 --with - mpi = 0 --download - f2cblaslapack = 1 --with - debugging
>> =
>> 1 [0]PETSC ERROR : #1 MatSetValues_SeqAIJ() line 485 in C :
>> \Apps\Projects\petsc - 3.6.3\src\mat\impls\aij\seq\aij.c [0]PETSC
>> ERROR : #2 MatSetValues() line 1173 in C : \Apps\Projects\petsc -
>> 3.6.3\src\mat\interface\matrix.c [0]PETSC ERROR : #3
>> MatDiagonalSet_Default() line 198 in C :
>> \Apps\Projects\petsc - 3.6.3\src\mat\utils\axpy.c [0]PETSC ERROR : #4
>> MatDiagonalSet_SeqAIJ() line 194 in C :
>> \Apps\Projects\petsc - 3.6.3\src\mat\impls\aij\seq\aij.c [0]PETSC
>> ERROR : #5 MatDiagonalSet() line 242 in C : \Apps\Projects\petsc -
>> 3.6.3\src\mat\utils\axpy.c [0]PETSC ERROR : #6
>> Tao_SSLS_FunctionGradient() line 66 in C :
>> \Apps\Projects\petsc -
>> 3.6.3\src\tao\complementarity\impls\ssls\ssls.c
>> [0]PETSC ERROR : #7 TaoLineSearchComputeObjectiveAndGradient() line
>> 941 in C
>> : \Apps\Projects\petsc -
>> 3.6.3\src\tao\linesearch\interface\taolinesearch.c
>> [0]PETSC ERROR : #8 TaoSolve_SSILS() line 64 in C :
>> \Apps\Projects\petsc -
>> 3.6.3\src\tao\complementarity\impls\ssls\ssils.c
>> [0]PETSC ERROR : #9 TaoSolve() line 204 in C : \Apps\Projects\petsc -
>> 3.6.3\src\tao\interface\taosolver.c
>> [0]PETSC ERROR : #10 run_complementarity() line 290 in c :
>> \apps\projects\complementarity.cpp
>>
>> I appreciate your help in this matter.
>>
>> Thanks
>>
>> Chris
>>
>> -----Original Message-----
>> From: Barry Smith [mailto:bsmith at mcs.anl.gov]
>> Sent: Monday, April 11, 2016 3:16 PM
>> To: Christopher Athaide <cnathaide at yahoo.com>
>> Cc: petsc-users at mcs.anl.gov
>> Subject: Re: [petsc-users] Problem setting Jacobian matrix in
>> complementarity problem
>>
>>
>>> On Apr 11, 2016, at 2:04 PM, Christopher Athaide
>>> <cnathaide at yahoo.com>
>> wrote:
>>>
>>> Hi,
>>>
>>> I am trying to solve a mixed complementarity problem using the TAO
>>> SSILS
>> Solver .
>>>
>>> C(x) = A x + d , l <= x <= u
>>>
>>> where A = n x n matrix
>>> D = n x 1 vector
>>>
>>> The Jacobian matrix is equal to the matrix A and is constant at each
>> iteration.
>>>
>>> I have setup the Jacobian update function
>>>
>>> PetscErrorCode FormJacobian(Tao, Vec X, Mat H, Mat, void *ptr) {
>>> AppCtx *data = (AppCtx *)ptr;
>>> PetscErrorCode ierr;
>>> PetscInt N;
>>> PetscInt zero = 0;
>>> PetscScalar *d, v;
>>> PetscBool assembled;
>>>
>>> PetscFunctionBegin;
>>> ierr = MatAssembled(H, &assembled); CHKERRQ(ierr);
>>> if (assembled) { ierr = MatZeroEntries(H);
>>> CHKERRQ(ierr); }
>>>
>>> /* Get pointers to vector data */
>>> VecGetSize(X, &N);
>>> for (auto i = zero; i != N; ++i)
>>> {
>>> for (auto j = zero; j != N; ++j)
>>> {
>>> ierr =
>> MatGetValue(data->A, i, j, &v);
>>> ierr = MatSetValue(H,
>>> i,
>> j, v, INSERT_VALUES);
>>> }
>>> }
>>>
>>> /* Assemble the matrix */
>>> ierr = MatSetOption(H, MAT_IGNORE_ZERO_ENTRIES,
>> PETSC_TRUE); CHKERRQ(ierr);
>>> ierr = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
>> CHKERRQ(ierr);
>>> ierr = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
>> CHKERRQ(ierr);
>>> PetscFunctionReturn(0);
>>>
>>> }
>>>
>>> In this case, after iterating over all the row and column indexes
>>> and
>> assembling the matrix the function seems to run correctly.
>>>
>>>
>>> But when I try to copy the matrix
>>>
>>> PetscErrorCode FormJacobian(Tao, Vec X, Mat H, Mat, void *ptr) {
>>>
>>> PetscFunctionReturn(0);
>>> AppCtx *data = (AppCtx *)ptr;
>>> PetscErrorCode ierr;
>>> ierr = MatAssembled(H, &assembled); CHKERRQ(ierr); if (assembled) {
>>> ierr = MatZeroEntries(H); CHKERRQ(ierr); } ierr =
>>> MatConvert(data->A, MATSAME, MatReuse::MAT_INITIAL_MATRIX, &H);
>>> CHKERRQ(ierr);
>>
>> You cannot do this because you are creating an entirely new H matrix
>> here in the subroutine that TAO won't know about.
>>
>> You should use MatCopy(data->A, H, SAME_NONZERO_PATTERN);
>>
>> But since H is constant why have a copy? Just set the values in H
>> initially and don't have a data->A matrix at all. Then your
>> FormJacobian() doesn't need to do anything since H is already set
>> with proper values and H remains the same through all iterations.
>>
>>
>>
>> Barry
>>
>>
>>>
>>>
>>>
>>> /* Assemble the matrix */
>>> ierr = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
>> CHKERRQ(ierr);
>>> ierr = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
>> CHKERRQ(ierr);
>>> PetscFunctionReturn(0); }
>>>
>>> I get the following error messages.
>>>
>>> Problem size: 259 rows
>>> [0]PETSC ERROR: --------------------- Error Message
>>> --------------------------------------------------------------
>>> [0]PETSC ERROR: Object is in wrong state [0]PETSC ERROR: Not for
>>> unassembled matrix [0]PETSC ERROR: See
>>> http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble
shooting.
>>> [0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015 [0]PETSC
>>> ERROR: Unknown Name on a arch-mswin-c-debug named LAGRANGE by
>>> Christopher Mon Apr 11 19:57:37 2016 [0]PETSC ERROR: Configure
>>> options --with-cc="win32fe cl" --with-fc=0 --with-mpi=0
>>> --download-f2cblaslapack=1 --with-debugging=1 [0]PETSC ERROR: #1
>>> MatMult() line 2210 in
>>> C:\Apps\Projects\petsc-3.6.3\src\mat\interface\matrix.c
>>> [0]PETSC ERROR: #2 MatDFischer() line 284 in C:\Apps\Projects
>>> \petsc-3.6.3\src\tao\util\tao_util.c
>>> [0]PETSC ERROR: #3 Tao_SSLS_FunctionGradient() line 64 in
>>> C:\Apps\Projects
>>> \petsc-3.6.3\src\tao\complementarity\impls\ssls\ssls.c
>>> [0]PETSC ERROR: #4 TaoLineSearchComputeObjectiveAndGradient() line
>>> 941 in C:\Apps\Projects
>>> \petsc-3.6.3\src\tao\linesearch\interface\taolinesearch.c
>>> [0]PETSC ERROR: #5 TaoSolve_SSILS() line 64 in C:\Apps\Projects
>>> \petsc-3.6.3\src\tao\complementarity\impls\ssls\ssils.c
>>> [0]PETSC ERROR: #6 TaoSolve() line 204 in C:\Apps\Projects\
>>> petsc-3.6.3\src\tao\interface\taosolver.c
>>> [0]PETSC ERROR: #7 run_complementarity() line 394 in
>>> c:\apps\projects \complementarity.cpp
>>>
>>> I am new to PETSC. Can anyone help me?
>>>
>>> Thanks
>>>
>>> Chris Athiade
>>
>
More information about the petsc-users
mailing list