[petsc-users] Problem with MatCreateSeqAIJWithArrays

[Chung-Kan Huang]

Hi,

I am writing a testing code that read CRS matrix M and rhs vector B from
two files and solve M*B=X.

I got error message indicating that I have segmentation fault but problem
disappeared after I removed MatDestroy(& M);

Another problem is that I am not sure if I should have
MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);

MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);

after

MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, nrow, ncol, rowptr, colptr,
mat, &M);

I attached my test code and hope someone can give me some pointers.

Thanks,

Kan


#include "petscksp.h"

#include "petscao.h"

#include "petsc.h"

#include "HYPRE.h"

#include "boost/timer/timer.hpp"

#include <string>

#include <iostream>

typedef boost::timer::cpu_timer Timer;

void P_Solve(const char * matrixFileName,

 const char * rhsFileName,

 const char * solFileName)

{

 Timer timer;

 PetscInitialize(PETSC_NULL,

 PETSC_NULL,

 PETSC_NULL,

 PETSC_NULL);

 Mat M; // initial matrix downloaded from disk

 Vec B; // initial RHS vector downloaded from disk

 Vec X; // initial guess=0/solution vector

 PC prec;

 KSP solver;

 FILE * fp;

 char temp[3];

 PetscInt nrow, ncol, nnz;

 PetscFOpen(PETSC_COMM_WORLD, matrixFileName, "r", &fp);

 fscanf(fp, "%s", temp); //#

 fscanf(fp, "%i %i %i", &nrow, &ncol, &nnz);

  PetscInt * rowptr = new PetscInt [nrow + 1];

 PetscInt * colptr = new PetscInt [nnz];

 PetscReal * mat = new PetscReal [nnz];

  fscanf(fp, "%s", temp); //#

 for (int i = 0; i < nrow + 1; i++) {

 fscanf(fp, "%i", &rowptr[i]);

 }

 fscanf(fp, "%s", temp); //#

 for (int i = 0; i < nnz; i++) {

 fscanf(fp, "%i", &colptr[i]);

 }

 fscanf(fp, "%s", temp); //#

 for (int i = 0; i < nnz; i++) {

 fscanf(fp, "%lf", &mat[i]);

 }

 fclose(fp);

  MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, nrow, ncol, rowptr, colptr,
mat, &M);

 MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);

 MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);

  // MatView(M, PETSC_VIEWER_STDOUT_WORLD);

 PetscFOpen(PETSC_COMM_WORLD, rhsFileName, "r", &fp);

  PetscReal * rhs = new PetscReal [nrow];

 PetscReal * x = new PetscReal [nrow];

  for (int i = 0; i < nrow; i++) {

 fscanf(fp, "%lf", &rhs[i]);

 }

 fclose(fp);

 VecCreateSeqWithArray(PETSC_COMM_WORLD, 1, nrow, rhs, &B);

 // VecView(B, PETSC_VIEWER_STDOUT_WORLD);

 MatGetVecs(M, &X, NULL);

 KSPCreate(PETSC_COMM_WORLD, &solver);

 KSPSetFromOptions(solver);

  KSPConvergedReason reason;

  PetscReal rtol = 0;

 PetscReal atol = 0;

 PetscReal dtol = 0;

 PetscInt it = 0;

 PetscInt its = 0;

 KSPGetTolerances( solver, &rtol, &atol, &dtol, &it);

 rtol = 1e-6;

 it = 500;

 KSPSetTolerances( solver, rtol, atol, dtol, it);

  KSPSetOperators(solver, M, M, SAME_NONZERO_PATTERN);

  KSPSetType(solver, KSPBCGS);

 KSPSetPCSide(solver, PC_RIGHT); // right preconditioning

 KSPSetInitialGuessNonzero(solver,PETSC_FALSE);

  KSPGetPC(solver, &prec);

 PCSetType(prec, PCILU);

 KSPSetFromOptions(solver);

 KSPSetUp(solver);

 PCSetUp(prec);

 PCFactorGetMatrix(prec, &M);

 timer.resume();

 KSPSolve(solver,B,X);

 timer.stop();

 KSPGetConvergedReason(solver,&reason);

 if (reason==KSP_DIVERGED_INDEFINITE_PC) {

 PetscPrintf(PETSC_COMM_WORLD,"\nDivergence because of indefinite
preconditioner;\n");

 PetscPrintf(PETSC_COMM_WORLD,"Run the executable again but with
'-pc_factor_shift_type POSITIVE_DEFINITE' option.\n");

 } else if (reason<0) {

 PetscPrintf(PETSC_COMM_WORLD,"\nOther kind of divergence: this should not
happen. %d\n",reason);

 }

 KSPGetIterationNumber(solver, &its);

 std::cout << "PETSc ITS = " << its << " "

 << "TIMES(SEC) = " << timer.elapsed().wall / 1.e9 << "\n";

 // VecView(X, PETSC_VIEWER_STDOUT_WORLD);

 VecGetArray(X, &x);

 PetscFOpen(PETSC_COMM_WORLD, solFileName, "w", &fp);

 for (int i = 0; i < nrow; i++) {

 fprintf(fp, "%.20e\n", x[i]);

 }

 VecRestoreArray(X, &x);

 fclose(fp);

 KSPDestroy(& solver);

 VecDestroy(& B);

 VecDestroy(& X);

 MatDestroy(& M);

 if (rowptr == NULL) delete [] rowptr;

 if (colptr == NULL) delete [] colptr;

 if (mat == NULL) delete [] mat;

 if (rhs == NULL) delete [] rhs;

   PetscFinalize();

}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130919/50446f01/attachment-0001.html>