# [petsc-users] Trying to solve Ax = b problem with Boundary condition using ksp

Matthew Knepley knepley at gmail.com
Thu Jul 12 06:22:31 CDT 2018

```On Thu, Jul 12, 2018 at 1:00 AM Fazlul Huq <huq2090 at gmail.com> wrote:

> Hi,
> I am trying to solve this problem with the following code:
>
> static char help[] = "Solves a tridiagonal linear system with KSP.\n\n";
>
> /*T
>    Concepts: KSP^solving a system of linear equations
>    Processors: 1
> 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 parallel example is ex23.c
> */
> #include <petscksp.h>
>
> int main(int argc,char **args)
> {
>   Vec            x, b, u;      /* approx solution, RHS, exact solution */
>   Mat            A;            /* linear system matrix */
>   KSP            ksp;          /* linear solver context */
>   PC             pc;           /* preconditioner context */
>   PetscReal      norm;         /* norm of solution error */
>   PetscErrorCode ierr;
>   PetscInt       i,n = 10,col[3],its;
>   PetscMPIInt    size;
>   // PetscScalar    zero = 0;
>   PetscScalar    hundredth = 0.01,value[3];
>   PetscScalar    one = 1.0;
>   PetscScalar    leftbc = 10.01;
>   PetscScalar    rightbc = 15.01;
>   PetscBool      nonzeroguess = PETSC_FALSE,changepcside = PETSC_FALSE;
>
>   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
>   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
>   if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor
> example only!");
>   ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
>   ierr =
> PetscOptionsGetBool(NULL,NULL,"-nonzero_guess",&nonzeroguess,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.
>   */
>   ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
>   ierr = PetscObjectSetName((PetscObject) x, "Solution");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);
>
>   i=1;
>     ierr = VecSetValues(b, 1, &i, &leftbc, INSERT_VALUES);CHKERRQ(ierr);
>   for (i=2; i<n-1; i++){
>     ierr = VecSetValues(b, n-2, &i, &hundredth,
> INSERT_VALUES);CHKERRQ(ierr);
>   }
>   i=n-1;
>     ierr = VecSetValues(b, 1, &i, &rightbc, INSERT_VALUES);CHKERRQ(ierr);
>     ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
>     ierr = VecAssemblyEnd(b);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.
>   */
>   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
>   ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
>   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
>   ierr = MatSetUp(A);CHKERRQ(ierr);
>
>   /*
>      Assemble matrix
>   */
>   value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
>   for (i=1; i<n-1; i++) {
>     col[0] = i-1; col[1] = i; col[2] = i+1;
>     ierr   = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
>   }
>   i    = n - 1; col[0] = n - 2; col[1] = n - 1;
>   ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
>   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);
>   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
>
>   /*
>      Set exact solution; then compute right-hand-side vector.
>   */
>
>
>   ierr = VecSet(u,one);CHKERRQ(ierr);
>   // ierr = MatMult(A,u,b);CHKERRQ(ierr);
>
>   i=0;
>     ierr = VecSetValues(b, 1, &i, &leftbc, INSERT_VALUES);CHKERRQ(ierr);
>   for (i=1; i<n-1; i++){
>     ierr = VecSetValues(b, n-2, &i, &hundredth,
> INSERT_VALUES);CHKERRQ(ierr);
>

What is the line above intended to do? You are telling it to take n-2
values, but only giving it
one value (hundreth) so it is indexing far outside the array.

Matt

>   }
>   i=n-1;
>     ierr = VecSetValues(b, 1, &i, &rightbc, INSERT_VALUES);CHKERRQ(ierr);
>     ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
>     ierr = VecAssemblyEnd(b);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);CHKERRQ(ierr);
>
>   /*
>      Test if you can change the KSP side and type after they have been
> previously set
>   */
>   ierr =
> PetscOptionsGetBool(NULL,NULL,"-change_pc_side",&changepcside,NULL);CHKERRQ(ierr);
>   if (changepcside) {
>     ierr = KSPSetUp(ksp);CHKERRQ(ierr);
>     ierr = PetscOptionsInsertString(NULL,"-ksp_norm_type unpreconditioned
> -ksp_pc_side right");CHKERRQ(ierr);
>     ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
>     ierr = KSPSetUp(ksp);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-5,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);
>
>   if (nonzeroguess) {
>     PetscScalar p = .5;
>     ierr = VecSet(x,p);CHKERRQ(ierr);
>     ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);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,-1.0,u);CHKERRQ(ierr);
>   ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
>   ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
>   ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations
> %D\n",(double)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_view).
>   */
>   ierr = PetscFinalize();
>   return ierr;
> }
>
> But I got the following error message:
>
> [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/documentation/faq.html#valgrind
> [0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS
> X 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: See http://www.mcs.anl.gov/petsc/documentation/faq.html
> for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.8.3, Dec, 09, 2017
> [0]PETSC ERROR: ./ex1 on a arch-linux2-c-opt named huq2090 by huq2090 Wed
> Jul 11 23:48:10 2018
> [0]PETSC ERROR: Configure options
> --with-ssl=0 --with-debugging=no --with-pic=1 --with-shared-libraries=1
> -F77=mpif77 -F90=mpif90 -CFLAGS="-fPIC -fopenmp" -CXXFLAGS="-fPIC -fopenmp"
> -FFLAGS="-fPIC -fopenmp" -FCFLAGS="-fPIC -fopenmp" -F90FLAGS="-fPIC
> -fopenmp" -F77FLAGS="-fPIC -fopenmp"
> [0]PETSC ERROR: #1 User provided function() line 0 in  unknown file
> application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
> [unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=59
> :
> system msg for write_line failure : Bad file descriptor
>
> And again, when I tried to solve the same code with another pc, I got
> different error message attached herewith as an attachment.
>
> Can you please suggest me about what to do.
>
> Thanks.
>
> Best,
> Huq
>
>
>
>
> --
>
> Fazlul Huq
> Graduate Research Assistant
> Department of Nuclear, Plasma & Radiological Engineering (NPRE)
> University of Illinois at Urbana-Champaign (UIUC)
> E-mail: huq2090 at gmail.com
>

--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their