[petsc-users] Problem setting Jacobian matrix in complementarity problem
Christopher Athaide
cnathaide at yahoo.com
Tue Apr 12 09:32:02 CDT 2016
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