[petsc-users] Reading a rectangular system from a file on disk and solving system in parallel

[Jed Brown]

On Wed, Jan 23, 2013 at 12:28 AM, Gaurish Telang <gaurish108 at gmail.com>wrote:

> Hi
>
> I want to solve a  least squares problem min || Ax - b|| in parallel with
> PETSc.
>
> I have the matrix A and the column vector b stored in files which are in
> PETSc' s binary format.
>
> So far, I am able to read both these files and solve the least squares
> system with PETSc's lsqr routine on a single processor.
>
> To solve the least squares system in parallel, can someone guide me on how
> to distribute the matrix A and vector b among several processors after
>  reading A and b from the corresponding files on disk?
>

Did you try just not VecSetType and MatSetType to force SEQ
implementations? It's better to call MatSetFromOptions() and
VecSetFromOptions() so you can change the type at run-time.

src/ksp/ksp/examples/tutorials/ex10.c automatically loads the matrix in
parallel and solves non-square systems (use -ksp_type lsqr).


>
>
> For reference, I have pasted my code below, which works well on a single
> processor.
>
> Thank you,
>
> Gaurish
>
>
>
>
> -----------------------------------------------------------------------------------------------------
>
> static char help[] = "--";
> #include <petscksp.h>
> #include <petscsys.h>
> #include <stdio.h>
> #undef __FUNCT__
> #define __FUNCT__ "main"
> int main(int argc,char **args)
> {
>   Vec            x, b, residue;      /* approx solution, RHS, residual */
>   Mat            A;            /* linear system matrix */
>   KSP            ksp;         /* linear solver context */
>   PC             pc;           /* preconditioner context */
>   PetscErrorCode ierr;
>   PetscInt       m,n   ;      /* # number of rows and columns of the
> matrix read in*/
>   PetscViewer    fd  ;
>   PetscInt       size, its;
>   PetscScalar    norm, tol=1.e-5;
>   PetscBool      flg;
>   char           fmat[PETSC_MAX_PATH_LEN];     /* input file names */
>   char           frhs[PETSC_MAX_PATH_LEN];     /* input file names */
>
>   PetscInitialize(&argc,&args,(char *)0,help);
>   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 =
>  PetscOptionsGetString(PETSC_NULL,"-fmat",fmat,PETSC_MAX_PATH_LEN,&flg);
>    ierr =
> PetscOptionsGetString(PETSC_NULL,"-frhs",frhs,PETSC_MAX_PATH_LEN,&flg);
>
>
>   /* Read in the matrix */
>   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, fmat, FILE_MODE_READ, &fd
> );  CHKERRQ(ierr);
>   ierr = MatCreate(PETSC_COMM_WORLD,&A);
>   ierr = MatSetType(A,MATSEQAIJ);
>   ierr = MatLoad(A,fd);
>   ierr = PetscViewerDestroy(&fd);
>
>
>   /*Read in the right hand side*/
>   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, frhs, FILE_MODE_READ, &fd
> );  CHKERRQ(ierr);
>   ierr = VecCreate(PETSC_COMM_WORLD,&b);
>   ierr = VecSetType(b,VECSEQ);
>   ierr = VecLoad(b,fd);
>   ierr = PetscViewerDestroy(&fd);
>
>
>   /* Get the matrix size. Used for setting the vector sizes*/
>   ierr = MatGetSize(A , &m, &n);
>   printf("The size of the matrix read in is %d x %d\n", m , n);
>
>   ierr = VecCreateSeq(PETSC_COMM_WORLD,n, &x);
>   ierr = VecCreateSeq(PETSC_COMM_WORLD,m, &residue);
>
>
>   /* Set the solver type at the command-line */
>   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
>   ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
>   ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
>   ierr = PCSetType(pc,PCNONE);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);
>
>   /* Initial guess for the krylov method set to zero*/
>   PetscScalar p = 0;
>   ierr = VecSet(x,p);CHKERRQ(ierr);
>   ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
>   /* 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);
>
>
>   /*View the numerical solution and the residue vector*/
>   printf("-------------------\nThe numerical solution is x = \n");
>   ierr = VecView(x,PETSC_VIEWER_STDOUT_SELF);
>
>
>   printf("-------------------\nThe residue vector Ax - b is = \n");
>   ierr = MatMult(A,x,residue);CHKERRQ(ierr);
>   ierr = VecAXPY(residue,-1.0,b);CHKERRQ(ierr);
>
>   ierr = VecView(residue,PETSC_VIEWER_STDOUT_SELF);
>
>   /* Clean up */
>   ierr = VecDestroy(&x);CHKERRQ(ierr);
>   ierr = VecDestroy(&residue);CHKERRQ(ierr);
>   ierr = VecDestroy(&b);CHKERRQ(ierr);
>   ierr = MatDestroy(&A);CHKERRQ(ierr);
>   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
>
>   ierr = PetscFinalize();
>   return 0;
> }
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130123/0e63795e/attachment.html>