[petsc-users] MatCreateSchurComplement
Anton Popov
popov at uni-mainz.de
Tue Oct 1 11:35:27 CDT 2013
Dear all,
I guess MatCreateSchurComplement is somehow broken in PETSc 3.4.2
Attached is the modified ex23 from src/ksp/ksp/examples/tutorials
I just copied system matrix three times and created a Schur complement
object.
However mpirun -n 4 ./ex23 gives: Unknown Mat type given: schurcomplement!
[0]PETSC ERROR:
------------------------------------------------------------------------
[0]PETSC ERROR: MatSetType() line 63 in
/Users/anton/LIB/petsc-3.4.2/src/mat/interface/matreg.c
[0]PETSC ERROR: MatCreateSchurComplement() line 180 in
/Users/anton/LIB/petsc-3.4.2/src/ksp/ksp/utils/schurm.c
[0]PETSC ERROR: main() line 121 in src/ksp/ksp/examples/tutorials/ex23.c
[1]PETSC ERROR: --------------------- Error Message
------------------------------------
[1]PETSC ERROR: Unknown type. Check for miss-spelling or missing
external package needed for type:
see http://www.mcs.anl.gov/petsc/documentation/installation.html#external!
[1]PETSC ERROR: Unknown Mat type given: schurcomplement!
[1]PETSC ERROR:
------------------------------------------------------------------------
[1]PETSC ERROR: Petsc Release Version 3.4.2, Jul, 02, 2013
[1]PETSC ERROR: See docs/changes/index.html for recent updates.
[1]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
[1]PETSC ERROR: See docs/index.html for manual pages.
[1]PETSC ERROR:
------------------------------------------------------------------------
[1]PETSC ERROR: ./ex23 on a arch-darwin-c-debug named
mac11-005.geo.uni-mainz.de by anton Tue Oct 1 18:30:50 2013
[1]PETSC ERROR: Libraries linked from /opt/petsc/petsc-3.4.2-int32-debug/lib
[1]PETSC ERROR: Configure run at Mon Sep 30 12:17:59 2013
[1]PETSC ERROR: Configure options
--prefix=/opt/petsc/petsc-3.4.2-int32-debug --download-f-blas-lapack=1
--with-debugging=1 --with-valgrind=1 --with-valgrind-dir=/opt/valgrind
--COPTFLAGS="-g -O0" --FOPTFLAGS="-g -O0" --CXXOPTFLAGS="-g -O0"
--with-shared-libraries=0 --with-cc=mpicc --with-cxx=mpicxx
--with-fc=mpif90 --download-ml=1 --download-hypre=1 --download-blacs=1
--download-scalapack=1 --download-metis=1 --download-parmetis=1
--download-mumps=1 --download-superlu_dist=1
Similar error we have in our code if we use MatCreateSchurComplement
function explicitely.
Any idea?
Thanks,
Anton
-------------- next part --------------
static char help[] = "Solves a tridiagonal linear system.\n\n";
/*T
Concepts: KSP^basic parallel example;
Processors: n
T*/
/*
Include "petscksp.h" so that we can use KSP solvers. Note that this file
automatically includes:
petscsys.h - base PETSc routines petscvec.h - vectors
petscmat.h - matrices
petscis.h - index sets petscksp.h - Krylov subspace methods
petscviewer.h - viewers petscpc.h - preconditioners
Note: The corresponding uniprocessor example is ex1.c
*/
#include <petscksp.h>
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
Vec x, b, u; /* approx solution, RHS, exact solution */
Mat A; /* linear system matrix */
Mat L, D, U, S;
KSP ksp; /* linear solver context */
PC pc; /* preconditioner context */
PetscReal norm,tol=1.e-11; /* norm of solution error */
PetscErrorCode ierr;
PetscInt i,n = 10,col[3],its,rstart,rend,nlocal;
PetscScalar neg_one = -1.0,one = 1.0,value[3];
PetscInitialize(&argc,&args,(char*)0,help);
ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Compute the matrix and right-hand-side vector that define
the linear system, Ax = b.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create vectors. Note that we form 1 vector from scratch and
then duplicate as needed. For this simple case let PETSc decide how
many elements of the vector are stored on each processor. The second
argument to VecSetSizes() below causes PETSc to decide.
*/
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
ierr = VecDuplicate(x,&u);CHKERRQ(ierr);
/* Identify the starting and ending mesh points on each
processor for the interior part of the mesh. We let PETSc decide
above. */
ierr = VecGetOwnershipRange(x,&rstart,&rend);CHKERRQ(ierr);
ierr = VecGetLocalSize(x,&nlocal);CHKERRQ(ierr);
/*
Create matrix. When using MatCreate(), the matrix format can
be specified at runtime.
Performance tuning note: For problems of substantial size,
preallocation of matrix memory is crucial for attaining good
performance. See the matrix chapter of the users manual for details.
We pass in nlocal as the "local" size of the matrix to force it
to have the same parallel layout as the vector created above.
*/
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,nlocal,nlocal,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
/*
Assemble matrix.
The linear system is distributed across the processors by
chunks of contiguous rows, which correspond to contiguous
sections of the mesh on which the problem is discretized.
For matrix assembly, each processor contributes entries for
the part that it owns locally.
*/
if (!rstart) {
rstart = 1;
i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
}
if (rend == n) {
rend = n-1;
i = n-1; col[0] = n-2; col[1] = n-1; value[0] = -1.0; value[1] = 2.0;
ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
}
/* Set entries corresponding to the mesh interior */
value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
for (i=rstart; i<rend; i++) {
col[0] = i-1; col[1] = i; col[2] = i+1;
ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
}
/* Assemble the matrix */
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
// duplicate matrix for test purposes
ierr = MatDuplicate(A, MAT_COPY_VALUES, &L); CHKERRQ(ierr);
ierr = MatDuplicate(A, MAT_COPY_VALUES, &D); CHKERRQ(ierr);
ierr = MatDuplicate(A, MAT_COPY_VALUES, &U); CHKERRQ(ierr);
// just check Schur complement creation
ierr = MatCreateSchurComplement(A, A, L, U, D, &S); CHKERRQ(ierr);
// destroy test matrices
ierr = MatDestroy(&L); CHKERRQ(ierr);
ierr = MatDestroy(&D); CHKERRQ(ierr);
ierr = MatDestroy(&U); CHKERRQ(ierr);
ierr = MatDestroy(&S); CHKERRQ(ierr);
/*
Set exact solution; then compute right-hand-side vector.
*/
ierr = VecSet(u,one);CHKERRQ(ierr);
ierr = MatMult(A,u,b);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the linear solver and set various options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create linear solver context
*/
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
/*
Set operators. Here the matrix that defines the linear system
also serves as the preconditioning matrix.
*/
ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
/*
Set linear solver defaults for this problem (optional).
- By extracting the KSP and PC contexts from the KSP context,
we can then directly call any KSP and PC routines to set
various options.
- The following four statements are optional; all of these
parameters could alternatively be specified at runtime via
KSPSetFromOptions();
*/
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
/*
Set runtime options, e.g.,
-ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
These options will override those specified above as long as
KSPSetFromOptions() is called _after_ any other customization
routines.
*/
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Solve linear system
*/
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
/*
View solver info; we could instead use the option -ksp_view to
print this info to the screen at the conclusion of KSPSolve().
*/
ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Check solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Check the error
*/
ierr = VecAXPY(x,neg_one,u);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
if (norm > tol) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G, Iterations %D\n",norm,its);CHKERRQ(ierr);
}
/*
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
*/
ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
/*
Always call PetscFinalize() before exiting a program. This routine
- finalizes the PETSc libraries as well as MPI
- provides summary and diagnostic information if certain runtime
options are chosen (e.g., -log_summary).
*/
ierr = PetscFinalize();
return 0;
}
More information about the petsc-users
mailing list