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

[Gaurish Telang]

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?


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/9b1f1371/attachment.html>