<div dir="ltr"><br clear="all"><span lang="EN"><div>Hi,</div><div> </div><div>I am writing a testing code that read CRS matrix M and rhs vector B from two files and solve M*B=X.</div><div> </div><div>I got error message indicating that I have segmentation fault but problem disappeared after I removed <span lang="EN">MatDestroy(& M);</span></div>
<div><span lang="EN"></span> </div><div><span lang="EN">Another problem is that I am not sure if I should have </span></div></span><div>MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);<p> MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);</p>
<p>after </p></div><div><span lang="EN"><p> MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, nrow, ncol, rowptr, colptr, mat, &M);</p></span></div><div> </div><div>I attached my test code and hope someone can give me some pointers.</div>
<div> </div><div>Thanks,</div><div> </div><div>Kan</div><div> </div><div><span lang="EN"><span lang="EN"><p>#include "petscksp.h"</p><p>
</p><p>#include "petscao.h"</p><p>
</p><p>#include "petsc.h"</p><p>
</p><p>#include "HYPRE.h"</p><p>
</p><p>#include "boost/timer/timer.hpp"</p><p>
</p><p>#include <string></p><p>
</p><p>#include <iostream></p><p>
</p><p>typedef boost::timer::cpu_timer Timer;</p><p>
</p><p>void P_Solve(const char * matrixFileName,</p><p>
</p><p> const char * rhsFileName,</p><p>
</p><p> const char * solFileName)</p><p>
</p><p>{</p><p>
</p><p> Timer timer;</p><p>
</p><p> PetscInitialize(PETSC_NULL,</p><p>
</p><p> PETSC_NULL,</p><p>
</p><p> PETSC_NULL,</p><p>
</p><p> PETSC_NULL);</p><p>
</p><p> Mat M; // initial matrix downloaded from disk</p><p>
</p><p> Vec B; // initial RHS vector downloaded from disk</p><p>
</p><p> Vec X; // initial guess=0/solution vector</p><p>
</p><p> PC prec;</p><p>
</p><p> KSP solver;</p><p>
</p><p> FILE * fp;</p><p>
</p><p> char temp[3];</p><p>
</p><p> PetscInt nrow, ncol, nnz;</p><p>
</p><p> PetscFOpen(PETSC_COMM_WORLD, matrixFileName, "r", &fp);</p><p>
</p><p> fscanf(fp, "%s", temp); //#</p><p>
</p><p> fscanf(fp, "%i %i %i", &nrow, &ncol, &nnz);</p><p>
</p><p> </p><p>
</p><p> PetscInt * rowptr = new PetscInt [nrow + 1];</p><p>
</p><p> PetscInt * colptr = new PetscInt [nnz]; </p><p>
</p><p> PetscReal * mat = new PetscReal [nnz];</p><p>
</p><p> </p><p>
</p><p> fscanf(fp, "%s", temp); //#</p><p>
</p><p> for (int i = 0; i < nrow + 1; i++) {</p><p>
</p><p> fscanf(fp, "%i", &rowptr[i]);</p><p>
</p><p> }</p><p>
</p><p> fscanf(fp, "%s", temp); //#</p><p>
</p><p> for (int i = 0; i < nnz; i++) {</p><p>
</p><p> fscanf(fp, "%i", &colptr[i]);</p><p>
</p><p> }</p><p>
</p><p> fscanf(fp, "%s", temp); //#</p><p>
</p><p> for (int i = 0; i < nnz; i++) {</p><p>
</p><p> fscanf(fp, "%lf", &mat[i]);</p><p>
</p><p> }</p><p>
</p><p> fclose(fp);</p><p>
</p><p> </p><p>
</p><p> MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, nrow, ncol, rowptr, colptr, mat, &M);</p><p>
</p><p> MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);</p><p>
</p><p> MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);</p><p>
</p><p> </p><p>
</p><p> // MatView(M, PETSC_VIEWER_STDOUT_WORLD);</p><p>
</p><p> PetscFOpen(PETSC_COMM_WORLD, rhsFileName, "r", &fp);</p><p>
</p><p> </p><p>
</p><p> PetscReal * rhs = new PetscReal [nrow];</p><p>
</p><p> PetscReal * x = new PetscReal [nrow];</p><p>
</p><p> </p><p>
</p><p> for (int i = 0; i < nrow; i++) {</p><p>
</p><p> fscanf(fp, "%lf", &rhs[i]);</p><p>
</p><p> }</p><p>
</p><p> fclose(fp); </p><p>
</p><p> VecCreateSeqWithArray(PETSC_COMM_WORLD, 1, nrow, rhs, &B);</p><p>
</p><p> // VecView(B, PETSC_VIEWER_STDOUT_WORLD);</p><p>
</p><p> MatGetVecs(M, &X, NULL);</p><p>
</p><p> KSPCreate(PETSC_COMM_WORLD, &solver);</p><p>
</p><p> KSPSetFromOptions(solver);</p><p>
</p><p> </p><p>
</p><p> KSPConvergedReason reason;</p><p>
</p><p> </p><p>
</p><p> PetscReal rtol = 0;</p><p>
</p><p> PetscReal atol = 0;</p><p>
</p><p> PetscReal dtol = 0;</p><p>
</p><p> PetscInt it = 0;</p><p>
</p><p> PetscInt its = 0;</p><p>
</p><p> KSPGetTolerances( solver, &rtol, &atol, &dtol, &it);</p><p>
</p><p> rtol = 1e-6;</p><p>
</p><p> it = 500;</p><p>
</p><p> KSPSetTolerances( solver, rtol, atol, dtol, it);</p><p>
</p><p> </p><p>
</p><p> KSPSetOperators(solver, M, M, SAME_NONZERO_PATTERN);</p><p>
</p><p> </p><p>
</p><p> KSPSetType(solver, KSPBCGS);</p><p>
</p><p> KSPSetPCSide(solver, PC_RIGHT); // right preconditioning</p><p>
</p><p> KSPSetInitialGuessNonzero(solver,PETSC_FALSE);</p><p>
</p><p> </p><p>
</p><p> KSPGetPC(solver, &prec);</p><p>
</p><p> PCSetType(prec, PCILU);</p><p>
</p><p> KSPSetFromOptions(solver);</p><p>
</p><p> KSPSetUp(solver); </p><p>
</p><p> PCSetUp(prec);</p><p>
</p><p> PCFactorGetMatrix(prec, &M);</p><p>
</p><p> timer.resume();</p><p>
</p><p> KSPSolve(solver,B,X);</p><p>
</p><p> timer.stop();</p><p>
</p><p> KSPGetConvergedReason(solver,&reason);</p><p>
</p><p> if (reason==KSP_DIVERGED_INDEFINITE_PC) {</p><p>
</p><p> PetscPrintf(PETSC_COMM_WORLD,"\nDivergence because of indefinite preconditioner;\n");</p><p>
</p><p> PetscPrintf(PETSC_COMM_WORLD,"Run the executable again but with '-pc_factor_shift_type POSITIVE_DEFINITE' option.\n");</p><p>
</p><p> } else if (reason<0) {</p><p>
</p><p> PetscPrintf(PETSC_COMM_WORLD,"\nOther kind of divergence: this should not happen. %d\n",reason);</p><p>
</p><p> }</p><p>
</p><p> KSPGetIterationNumber(solver, &its);</p><p>
</p><p> std::cout << "PETSc ITS = " << its << " "</p><p>
</p><p> << "TIMES(SEC) = " << timer.elapsed().wall / 1.e9 << "\n";</p><p>
</p><p> // VecView(X, PETSC_VIEWER_STDOUT_WORLD);</p><p>
</p><p> VecGetArray(X, &x);</p><p>
</p><p> PetscFOpen(PETSC_COMM_WORLD, solFileName, "w", &fp);</p><p>
</p><p> for (int i = 0; i < nrow; i++) {</p><p>
</p><p> fprintf(fp, "%.20e\n", x[i]);</p><p>
</p><p> }</p><p>
</p><p> VecRestoreArray(X, &x);</p><p>
</p><p> fclose(fp);</p><p>
</p><p> KSPDestroy(& solver); </p><p>
</p><p> VecDestroy(& B);</p><p>
</p><p> VecDestroy(& X);</p><p>
</p><p> MatDestroy(& M);</p><p>
</p><p> if (rowptr == NULL) delete [] rowptr;</p><p>
</p><p> if (colptr == NULL) delete [] colptr;</p><p>
</p><p> if (mat == NULL) delete [] mat;</p><p>
</p><p> if (rhs == NULL) delete [] rhs;</p><p>
</p><p> </p><p>
</p><p> </p><p>
</p><p> PetscFinalize(); </p><p>
</p><p>}</p><p>
</p><p> </p><p>
</p></span><p></p></span><p></p></div><p> </p>
</div>