Manfred :<div>Thanks for sending the modified code. I can repeat the error in 'MatFactorNumeric_MUMPS'.</div><div>Running it with '-mat_type aij' works well, likely there is a bug in our code.</div><div>I'll work on it and get back to you about the fix.</div>
<div><br></div><div>Meanwhile, you may use aij format instead of sbaij. </div><div><br></div><div>Hong</div><div><br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<u></u>
<div bgcolor="#ffffff" text="#000000">
Jed Brown schrieb:
<div><div class="h5"><blockquote type="cite">
<div class="gmail_quote">On Fri, Mar 2, 2012 at 04:18, Manfred Gratt <span dir="ltr"><<a href="mailto:manfred.gratt@uibk.ac.at" target="_blank">manfred.gratt@uibk.ac.at</a>></span>
wrote:<br>
<blockquote class="gmail_quote" style="border-left:1px solid rgb(204,204,204);margin:0pt 0pt 0pt 0.8ex;padding-left:1ex">thank
you for your quick response. The petsc-dev did not solve the
segfault. I checked what the difference from ex5 to my code was and
after I changed my code to ex5 the segfault was gone when I used
MatZeroEntries instead of destroying and creating a new Matrix for the
second solve. This seems to be logical but, what I do not understand
is why it works with destroying and creating a new Matrix on more than
one processor fine?</blockquote>
<div><br>
</div>
<div>I can't understand from this description exactly what works and
what doesn't work for you. There shouldn't be a SEGV for any
"reasonable" sequence of calls you make, so we should clarify the
confusion and either fix the bug or make a better/earlier error message.</div>
</div>
<br>
<div>Is there a way you can modify ex5 (or another example, or send
your own code) to show the problem?</div>
</blockquote></div></div>
I tried to change the ex5 to bring the same error but I am not able to
get the same error. I still get an error although it is an error I
don't understand. <br>
I changed ex5 to calculate in second solve the same matrix as in the
first solve. As written in definition of SAME_NONZERO_PATTERN it should
work. <br>
When before the second solve the matrix C is set to zero with
MatZeroEntries it works well, but when I instead destroy the matrix C
and create a new C it does not work with the error that it is an
Numerically singular matrix. Should that happen? It is the same matrix
as in the first solve. So shouldnt it happen there?<br>
When I use instead DIFFERENT_NONZERO_PATTERN it also works well in the
second solve.<br>
<br>
I modified the ex5 so that you can easily call the changes with
-mat_zeroent to call MatZeroEntries instead of destroy and create a new
matrix.<br>
The second option -mat_different calls DIFFERENT_NONZERO_PATTERN
instead of SAME_NONZERO_PATTERN.<br>
<br>
Without these options it throws an error when you call it with<div class="im"><br>
./ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky
-pc_factor_mat_solver_package mumps<br>
<br></div>
Thanks<br>
Manfred<br>
<br>
newex5.c:<br>
<br>
static char help[] = "Solves two linear systems in parallel with KSP.
The code\n\<br>
illustrates repeated solution of linear systems with the same
preconditioner\n\<br>
method but different matrices (having the same nonzero structure). The
code\n\<br>
also uses multiple profiling stages. Input arguments are\n\<br>
-m <size> : problem size\n\<br>
-mat_nonsym : use nonsymmetric matrix (default is symmetric)\n\n";<br>
<br>
/*T<br>
Concepts: KSP^repeatedly solving linear systems;<br>
Concepts: PetscLog^profiling multiple stages of code;<br>
Processors: n<br>
T*/<br>
<br>
/* <br>
Include "petscksp.h" so that we can use KSP solvers. Note that this
file<br>
automatically includes:<br>
petscsys.h - base PETSc routines petscvec.h - vectors<br>
petscmat.h - matrices<br>
petscis.h - index sets petscksp.h - Krylov subspace
methods<br>
petscviewer.h - viewers petscpc.h - preconditioners<br>
*/<br>
#include <petscksp.h><br>
<br>
#undef __FUNCT__<br>
#define __FUNCT__ "main"<br>
int main(int argc,char **args)<br>
{<br>
KSP ksp; /* linear solver context */<br>
Mat C; /* matrix */<br>
Vec x,u,b; /* approx solution, RHS, exact
solution */<br>
PetscReal norm; /* norm of solution error */<br>
PetscScalar v,none = -1.0;<br>
PetscInt Ii,J,ldim,low,high,iglobal,Istart,Iend;<br>
PetscErrorCode ierr;<br>
PetscInt i,j,m = 3,n = 2,its;<br>
PetscMPIInt size,rank;<br>
PetscBool mat_nonsymmetric = PETSC_FALSE;<br>
PetscBool matzeroent = PETSC_FALSE, different_pattern =
PETSC_FALSE;<br>
#if defined (PETSC_USE_LOG)<br>
PetscLogStage stages[2];<br>
#endif<br>
<br>
PetscInitialize(&argc,&args,(char *)0,help);<br>
ierr =
PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);CHKERRQ(ierr);<br>
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);<br>
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);<br>
n = 2*size;<br>
<br>
/*<br>
Set flag if we are doing a nonsymmetric problem; the default is
symmetric.<br>
*/<br>
ierr =
PetscOptionsGetBool(PETSC_NULL,"-mat_nonsym",&mat_nonsymmetric,PETSC_NULL);CHKERRQ(ierr);<br>
ierr =
PetscOptionsGetBool(PETSC_NULL,"-mat_zeroent",&matzeroent,PETSC_NULL);CHKERRQ(ierr);<br>
ierr =
PetscOptionsGetBool(PETSC_NULL,"-mat_different",&different_pattern,PETSC_NULL);CHKERRQ(ierr);<br>
/*<br>
Register two stages for separate profiling of the two linear
solves.<br>
Use the runtime option -log_summary for a printout of performance<br>
statistics at the program's conlusion.<br>
*/<br>
ierr = PetscLogStageRegister("Original
Solve",&stages[0]);CHKERRQ(ierr);<br>
ierr = PetscLogStageRegister("Second
Solve",&stages[1]);CHKERRQ(ierr);<br>
<br>
/* -------------- Stage 0: Solve Original System
---------------------- */<br>
/* <br>
Indicate to PETSc profiling that we're beginning the first stage<br>
*/<br>
ierr = PetscLogStagePush(stages[0]);CHKERRQ(ierr);<br>
<br>
/* <br>
Create parallel matrix, specifying only its global dimensions.<br>
When using MatCreate(), the matrix format can be specified at<br>
runtime. Also, the parallel partitioning of the matrix is<br>
determined by PETSc at runtime.<br>
*/<br>
ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);<br>
ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);<br>
ierr = MatSetFromOptions(C);CHKERRQ(ierr);<br>
ierr = MatSetUp(C);CHKERRQ(ierr);<br>
<br>
/* <br>
Currently, all PETSc parallel matrix formats are partitioned by<br>
contiguous chunks of rows across the processors. Determine which<br>
rows of the matrix are locally owned. <br>
*/<br>
ierr = MatGetOwnershipRange(C,&Istart,&Iend);CHKERRQ(ierr);<br>
<br>
/* <br>
Set matrix entries matrix in parallel.<br>
- Each processor needs to insert only elements that it owns<br>
locally (but any non-local elements will be sent to the<br>
appropriate processor during matrix assembly). <br>
- Always specify global row and columns of matrix entries.<br>
*/<br>
for (Ii=Istart; Ii<Iend; Ii++) { <br>
v = -1.0; i = Ii/n; j = Ii - i*n; <br>
if (i>0) {J = Ii - n; ierr =
MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}<br>
if (i<m-1) {J = Ii + n; ierr =
MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}<br>
if (j>0) {J = Ii - 1; ierr =
MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}<br>
if (j<n-1) {J = Ii + 1; ierr =
MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}<br>
v = 4.0; ierr =
MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);<br>
}<br>
<br>
/*<br>
Make the matrix nonsymmetric if desired<br>
*/<br>
if (mat_nonsymmetric) {<br>
for (Ii=Istart; Ii<Iend; Ii++) { <br>
v = -1.5; i = Ii/n;<br>
if (i>1) {J = Ii-n-1; ierr =
MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}<br>
}<br>
} else {<br>
ierr = MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);<br>
ierr =
MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr);<br>
}<br>
<br>
/* <br>
Assemble matrix, using the 2-step process:<br>
MatAssemblyBegin(), MatAssemblyEnd()<br>
Computations can be done while messages are in transition<br>
by placing code between these two statements.<br>
*/<br>
ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br>
ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br>
<br>
/* <br>
Create parallel vectors.<br>
- When using VecSetSizes(), we specify only the vector's global<br>
dimension; the parallel partitioning is determined at runtime. <br>
- Note: We form 1 vector from scratch and then duplicate as
needed.<br>
*/<br>
ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);<br>
ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);<br>
ierr = VecSetFromOptions(u);CHKERRQ(ierr);<br>
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);<br>
ierr = VecDuplicate(b,&x);CHKERRQ(ierr);<br>
<br>
/* <br>
Currently, all parallel PETSc vectors are partitioned by<br>
contiguous chunks across the processors. Determine which<br>
range of entries are locally owned.<br>
*/<br>
ierr = VecGetOwnershipRange(x,&low,&high);CHKERRQ(ierr);<br>
<br>
/*<br>
Set elements within the exact solution vector in parallel.<br>
- Each processor needs to insert only elements that it owns<br>
locally (but any non-local entries will be sent to the<br>
appropriate processor during vector assembly).<br>
- Always specify global locations of vector entries.<br>
*/<br>
ierr = VecGetLocalSize(x,&ldim);CHKERRQ(ierr);<br>
for (i=0; i<ldim; i++) {<br>
iglobal = i + low;<br>
v = (PetscScalar)(i + 100*rank);<br>
ierr =
VecSetValues(u,1,&iglobal,&v,INSERT_VALUES);CHKERRQ(ierr);<br>
}<br>
<br>
/* <br>
Assemble vector, using the 2-step process:<br>
VecAssemblyBegin(), VecAssemblyEnd()<br>
Computations can be done while messages are in transition,<br>
by placing code between these two statements.<br>
*/<br>
ierr = VecAssemblyBegin(u);CHKERRQ(ierr);<br>
ierr = VecAssemblyEnd(u);CHKERRQ(ierr);<br>
<br>
/* <br>
Compute right-hand-side vector<br>
*/<br>
ierr = MatMult(C,u,b);CHKERRQ(ierr);<br>
<br>
/* <br>
Create linear solver context<br>
*/<br>
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);<br>
<br>
/* <br>
Set operators. Here the matrix that defines the linear system<br>
also serves as the preconditioning matrix.<br>
*/<br>
if(different_pattern) {<br>
ierr =
KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); <br>
}else {<br>
ierr = KSPSetOperators(ksp,C,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr);<br>
}<br>
/* <br>
Set runtime options (e.g., -ksp_type <type> -pc_type
<type>)<br>
*/<br>
<br>
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);<br>
<br>
/* <br>
Solve linear system. Here we explicitly call KSPSetUp() for more<br>
detailed performance monitoring of certain preconditioners, such<br>
as ICC and ILU. This call is optional, as KSPSetUp() will<br>
automatically be called within KSPSolve() if it hasn't been<br>
called already.<br>
*/<br>
ierr = KSPSetUp(ksp);CHKERRQ(ierr);<br>
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);<br>
<br>
/* <br>
Check the error<br>
*/<br>
// ierr = VecAXPY(x,none,u);CHKERRQ(ierr);<br>
// ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);<br>
// ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);<br>
// if (norm > 1.e-13){<br>
// ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G,
Iterations %D\n",norm,its);CHKERRQ(ierr);<br>
// }<br>
<br>
/* -------------- Stage 1: Solve Second System ----------------------
*/<br>
/* <br>
Solve another linear system with the same method. We reuse the KSP<br>
context, matrix and vector data structures, and hence save the<br>
overhead of creating new ones.<br>
<br>
Indicate to PETSc profiling that we're concluding the first<br>
stage with PetscLogStagePop(), and beginning the second stage with<br>
PetscLogStagePush().<br>
*/<br>
ierr = PetscLogStagePop();CHKERRQ(ierr);<br>
ierr = PetscLogStagePush(stages[1]);CHKERRQ(ierr);<br>
<br>
/* <br>
Initialize all matrix entries to zero. MatZeroEntries() retains
the<br>
nonzero structure of the matrix for sparse formats.<br>
*/<br>
if(matzeroent) {<br>
ierr = MatZeroEntries(C);CHKERRQ(ierr);<br>
} else {<br>
ierr = MatDestroy(&C);<br>
<br>
ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);<br>
ierr =
MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);<br>
ierr = MatSetFromOptions(C);CHKERRQ(ierr);<br>
ierr = MatSetUp(C);CHKERRQ(ierr);<br>
}<br>
/* <br>
Assemble matrix again. Note that we retain the same matrix data<br>
structure and the same nonzero pattern; we just change the values<br>
of the matrix entries.<br>
*/<br>
for (Ii=Istart; Ii<Iend; Ii++) { <br>
v = -1.0; i = Ii/n; j = Ii - i*n; <br>
if (i>0) {J = Ii - n; ierr =
MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}<br>
if (i<m-1) {J = Ii + n; ierr =
MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}<br>
if (j>0) {J = Ii - 1; ierr =
MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}<br>
if (j<n-1) {J = Ii + 1; ierr =
MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}<br>
v = 4.0; ierr =
MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);<br>
}<br>
if (mat_nonsymmetric) {<br>
for (Ii=Istart; Ii<Iend; Ii++) { <br>
v = -1.5; i = Ii/n;<br>
if (i>1) {J = Ii-n-1; ierr =
MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}<br>
}<br>
} else {<br>
ierr = MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);<br>
ierr =
MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr);<br>
}<br>
ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br>
ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); <br>
<br>
/* <br>
Compute another right-hand-side vector<br>
*/<br>
ierr = MatMult(C,u,b);CHKERRQ(ierr);<br>
<br>
/* <br>
Set operators. Here the matrix that defines the linear system<br>
also serves as the preconditioning matrix.<br>
- The flag SAME_NONZERO_PATTERN indicates that the<br>
preconditioning matrix has identical nonzero structure<br>
as during the last linear solve (although the values of<br>
the entries have changed). Thus, we can save some<br>
work in setting up the preconditioner (e.g., no need to<br>
redo symbolic factorization for ILU/ICC preconditioners).<br>
- If the nonzero structure of the matrix is different during<br>
the second linear solve, then the flag DIFFERENT_NONZERO_PATTERN<br>
must be used instead. If you are unsure whether the<br>
matrix structure has changed or not, use the flag<br>
DIFFERENT_NONZERO_PATTERN.<br>
- Caution: If you specify SAME_NONZERO_PATTERN, PETSc<br>
believes your assertion and does not check the structure<br>
of the matrix. If you erroneously claim that the structure<br>
is the same when it actually is not, the new preconditioner<br>
will not function correctly. Thus, use this optimization<br>
feature with caution!<br>
*/<br>
<br>
if(different_pattern) {<br>
ierr =
KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); <br>
}else {<br>
ierr = KSPSetOperators(ksp,C,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr);<br>
}<br>
<br>
/* <br>
Solve linear system<br>
*/<br>
ierr = KSPSetUp(ksp);CHKERRQ(ierr);<br>
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); <br>
<br>
/* <br>
Check the error<br>
*/<br>
ierr = VecAXPY(x,none,u);CHKERRQ(ierr);<br>
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);<br>
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);<br>
if (norm > 1.e-4){<br>
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G, Iterations
%D\n",norm,its);CHKERRQ(ierr);<br>
}<br>
<br>
/* <br>
Free work space. All PETSc objects should be destroyed when they<br>
are no longer needed.<br>
*/<br>
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);<br>
ierr = VecDestroy(&u);CHKERRQ(ierr);<br>
ierr = VecDestroy(&x);CHKERRQ(ierr);<br>
ierr = VecDestroy(&b);CHKERRQ(ierr);<br>
ierr = MatDestroy(&C);CHKERRQ(ierr);<br>
<br>
/*<br>
Indicate to PETSc profiling that we're concluding the second stage
<br>
*/<br>
ierr = PetscLogStagePop();CHKERRQ(ierr);<br>
<br>
ierr = PetscFinalize();<br>
return 0;<br>
}<br>
<br>
<br>
<br>
<br>
<br>
</div>
</blockquote></div><br></div>