Direct LU solver

Thu Feb 28 13:07:36 CST 2008

Hi,

I need to do direct LU solves (repeatedly, with the same matrix) in one of
my MPI applications. I am having trouble implementing the solver. To
identify the problem, I wrote a short toy  code to run with 2 processes. I
can run it with either the spooles parallel matrix or the superlu_dist
matrix. I am using C++ and complex matrices.

Here is the code listing:

#include <iostream>
#include <complex>
#include <stdlib.h>
#include "petsc.h"
#include "petscmat.h"
#include "petscvec.h"
#include "petscksp.h"

using namespace std;

int main( int argc, char *argv[] ) {

  int rank, size;
  Mat A;
  PetscErrorCode ierr;
  PetscInt loc;
  PetscScalar val;
  Vec x, y;
  KSP solver;
  PC prec;
  MPI_Comm comm;

  // Number of non-zeros in each row
  int d_nnz=1, o_nnz=1;

  ierr=PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
  // Initialization (including MPI)

  comm=PETSC_COMM_WORLD;

  ierr=MPI_Comm_size(PETSC_COMM_WORLD,&size); CHKERRQ(ierr);
  ierr=MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);

  // Assemble matrix A

  if(rank==0) {
    ierr=MatCreateMPIAIJ(comm,1,1,2,2,0,&d_nnz,0,&o_nnz,&A); CHKERRQ(ierr);
    val=complex<double>(1.0,0.0);
    ierr=MatSetValue(A,0,0,val,INSERT_VALUES);CHKERRQ(ierr);
    val=complex<double>(0.0,1.0);
    ierr=MatSetValue(A,0,1,val,INSERT_VALUES);CHKERRQ(ierr);
  }
  else {
    ierr=MatCreateMPIAIJ(PETSC_COMM_WORLD,1,1,2,2,0,&d_nnz,0,&o_nnz,&A);
CHKERRQ(ierr);
    val=complex<double>(1.0,1.0);
    ierr=MatSetValue(A,1,1,val,INSERT_VALUES);CHKERRQ(ierr);
    val=complex<double>(0.0,-1.0);
    ierr=MatSetValue(A,1,0,val,INSERT_VALUES);CHKERRQ(ierr);
  }

  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

  cout << "============  Mat A  ==================" << endl;
  ierr=MatView(A,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  cout << "======================================" << endl;

  // For spooles
  //ierr=MatConvert(A,MATMPIAIJSPOOLES,MAT_REUSE_MATRIX,&A); CHKERRQ(ierr);

  // For superlu_dist
  ierr=MatConvert(A,MATSUPERLU_DIST,MAT_REUSE_MATRIX,&A); CHKERRQ(ierr);

  // Direct LU solver
  ierr=KSPCreate(comm,&solver); CHKERRQ(ierr);
  ierr=KSPSetType(solver,KSPPREONLY); CHKERRQ(ierr);
  ierr=KSPSetOperators(solver,A,A,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
  ierr=KSPGetPC(solver,&prec); CHKERRQ(ierr);
  ierr=PCSetType(prec,PCLU); CHKERRQ(ierr);
  ierr=KSPSetFromOptions(solver); CHKERRQ(ierr);

  //============  Vector assembly ========================

  if(rank==0) {
    ierr=VecCreateMPI(PETSC_COMM_WORLD,1,2,&x); CHKERRQ(ierr);
    val=complex<double>(1.0,0.0);
    loc=0;
    ierr=VecSetValues(x,1,&loc,&val,ADD_VALUES);CHKERRQ(ierr);
  }
  else {
    ierr=VecCreateMPI(PETSC_COMM_WORLD,1,2,&x); CHKERRQ(ierr);
    val=complex<double>(-1.0,0.0);
    loc=1;
    ierr=VecSetValues(x,1,&loc,&val,ADD_VALUES);CHKERRQ(ierr);
  }

  ierr=VecAssemblyBegin(x); CHKERRQ(ierr);
  ierr=VecAssemblyEnd(x); CHKERRQ(ierr);

  cout << "============== Vec x ==================" << endl;
  ierr=VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  cout << "======================================" << endl;

  VecDuplicate(x,&y);  // Duplicate the matrix storage

  // Solve the matrix equation
  ierr=KSPSolve(solver,x,y); CHKERRQ(ierr);

  cout << "============== Vec y =================" << endl;
  ierr=VecView(y,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  cout << "======================================" << endl;

  // Destructors
  ierr=KSPDestroy(solver); CHKERRQ(ierr);
  ierr=VecDestroy(x); CHKERRQ(ierr);
  ierr=VecDestroy(y); CHKERRQ(ierr);
  ierr=MatDestroy(A); CHKERRQ(ierr);

 // Finalize
  ierr=PetscFinalize(); CHKERRQ(ierr);

  return 0;

}

When I run the program with spooles, I get the following output.

============  Mat A  ==================
============  Mat A  ==================
======================================
row 0: (0, 1)  (1, 0 + 1 i)
row 1: (0, 0 - 1 i) (1, 1 + 1 i)
======================================
============== Vec x ==================
============== Vec x ==================
Process [0]
1
======================================
Process [1]
-1
======================================

 fatal error in InpMtx_MPI_split()
 firsttag = 0, tagbound = -1

 fatal error in InpMtx_MPI_split()
 firsttag = 0, tagbound = -1
-----------------------------------------------------------------------------
One of the processes started by mpirun has exited with a nonzero exit
code.  This typically indicates that the process finished in error.
If your process did not finish in error, be sure to include a "return
0" or "exit(0)" in your C code before exiting the application.

PID 22881 failed on node n0 (127.0.0.1) with exit status 255.
-----------------------------------------------------------------------------
mpirun failed with exit status 255

When run in the debugger, there is no stack trace. The error is

Program exited with code 0377

With superlu_dist, the output is

============  Mat A  ==================
============  Mat A  ==================
row 0: (0, 1)  (1, 0 + 1 i)
row 1: (0, 0 - 1 i) (1, 1 + 1 i)
======================================
======================================
============== Vec x ==================
============== Vec x ==================
Process [0]
1
Process [1]
-1
======================================
======================================
[0]PETSC ERROR:
------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation,
probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see
http://www.mcs.anl.gov/petsc/petsc-as/documentation/troubleshooting.html#Signal[0]PETSC
 ERROR: or try http://valgrind.org on linux or man libgmalloc on Apple to
find memory corruption errors
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and
run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: --------------------- Error Message
------------------------------------
[0]PETSC ERROR: Signal received!
[0]PETSC ERROR:
------------------------------------------------------------------------
[0]PETSC ERROR: Petsc Release Version 2.3.3, Patch 8, Fri Nov 16 17:03:40
CST 2007 HG revision: 414581156e67e55c761739b0deb119f7590d0f4b
[0]PETSC ERROR: See docs/changes/index.html for recent updates.
[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
[0]PETSC ERROR: See docs/index.html for manual pages.
[0]PETSC ERROR:
------------------------------------------------------------------------
[0]PETSC ERROR: lu on a linux-gnu named tabla by amit Thu Feb 28 13:48:54
2008
[0]PETSC ERROR: Libraries linked from
/home/amit/programs/ParEM/petsc-2.3.3-p8/lib/linux-gnu-c-debug
[0]PETSC ERROR: Configure run at Thu Feb 28 12:19:39 2008
[0]PETSC ERROR: Configure options --with-scalar-type=complex
--with-debugging=no --with-fortran-kernels=generic --with-clanguage=cxx
--with-metis=1 --download-metis=1 --with-parmetis=1 --download-parmetis=1
--with-superlu_dist=1 --download-superlu_dist=1 --with-spooles=1
--with-spooles-dir=/home/amit/programs/ParEM/spooles-2.2 COPTFLAGS="-O3
-march=p4 -mtune=p4 -ffast-math -malign-double -funroll-loops -pipe
-fomit-frame-pointer -finline-functions -msse2" CXXOPTFLAGS="-O3 -march=p4
-mtune=p4 -ffast-math -malign-double -funroll-loops -pipe
-fomit-frame-pointer -finline-functions -msse2" FOPTS="-O3 -qarch=p4
-qtune=p4" --with-shared=0
[0]PETSC ERROR:
------------------------------------------------------------------------
[0]PETSC ERROR: User provided function() line 0 in unknown directory
unknown file
-----------------------------------------------------------------------------
One of the processes started by mpirun has exited with a nonzero exit
code.  This typically indicates that the process finished in error.
If your process did not finish in error, be sure to include a "return
0" or "exit(0)" in your C code before exiting the application.

PID 22998 failed on node n0 (127.0.0.1) with exit status 1.
-----------------------------------------------------------------------------
mpirun failed with exit status 1

The debugger tracks the segmentation violation to

main -> KSPSolve -> KSPSetUp -> PCSetUp -> PCSetUp_LU -> MatLUFactorNumeric
-> MatLUFactorNumeric_SUperLU_DIST -> pzgssvx

Could  someone kindly point out what I am missing ?

Thanks

Rgds,
Amit