[petsc-users] Problem setting Jacobian matrix in complementarity problem
Barry Smith
bsmith at mcs.anl.gov
Tue Apr 12 13:07:50 CDT 2016
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