Hi <div><br></div><div>I want to solve a  least squares problem min || Ax - b|| in parallel with PETSc. </div><div><br></div><div>I have the matrix A and the column vector b stored in files which are in PETSc' s binary format.</div>
<div><br></div><div>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. </div><div><br></div><div>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?</div>
<div><br></div><div><br></div><div>For reference, I have pasted my code below, which works well on a single processor.</div><div><br></div><div>Thank you,</div><div><br></div><div>Gaurish</div><div><br></div><div><br></div>
<div><br></div><div>-----------------------------------------------------------------------------------------------------</div><div><br></div><div><div>static char help[] = "--";</div><div>#include <petscksp.h></div>
<div>#include <petscsys.h></div><div>#include <stdio.h></div><div>#undef __FUNCT__</div><div>#define __FUNCT__ "main"</div><div>int main(int argc,char **args)</div><div>{</div><div>  Vec            x, b, residue;      /* approx solution, RHS, residual */</div>
<div>  Mat            A;            /* linear system matrix */</div><div>  KSP            ksp;         /* linear solver context */</div><div>  PC             pc;           /* preconditioner context */</div><div>  PetscErrorCode ierr;</div>
<div>  PetscInt       m,n   ;      /* # number of rows and columns of the matrix read in*/</div><div>  PetscViewer    fd  ;</div><div>  PetscInt       size, its;</div><div>  PetscScalar    norm, tol=1.e-5;</div><div>  PetscBool      flg;</div>
<div>  char           fmat[PETSC_MAX_PATH_LEN];     /* input file names */</div><div>  char           frhs[PETSC_MAX_PATH_LEN];     /* input file names */</div><div><br></div><div>  PetscInitialize(&argc,&args,(char *)0,help);</div>
<div>  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);</div><div>  if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only!");</div><div> </div><div><br></div><div>   ierr =  PetscOptionsGetString(PETSC_NULL,"-fmat",fmat,PETSC_MAX_PATH_LEN,&flg);</div>
<div>   ierr = PetscOptionsGetString(PETSC_NULL,"-frhs",frhs,PETSC_MAX_PATH_LEN,&flg);</div><div><br></div><div><br></div><div>  /* Read in the matrix */</div><div>  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, fmat, FILE_MODE_READ, &fd );  CHKERRQ(ierr);</div>
<div>  ierr = MatCreate(PETSC_COMM_WORLD,&A);</div><div>  ierr = MatSetType(A,MATSEQAIJ);</div><div>  ierr = MatLoad(A,fd);</div><div>  ierr = PetscViewerDestroy(&fd);</div><div><br></div><div><br></div><div>  /*Read in the right hand side*/</div>
<div>  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, frhs, FILE_MODE_READ, &fd );  CHKERRQ(ierr);</div><div>  ierr = VecCreate(PETSC_COMM_WORLD,&b);</div><div>  ierr = VecSetType(b,VECSEQ);</div><div>  ierr = VecLoad(b,fd);</div>
<div>  ierr = PetscViewerDestroy(&fd);</div><div><br></div><div><br></div><div>  /* Get the matrix size. Used for setting the vector sizes*/</div><div>  ierr = MatGetSize(A , &m, &n);</div><div>  printf("The size of the matrix read in is %d x %d\n", m , n);</div>
<div>    </div><div>  ierr = VecCreateSeq(PETSC_COMM_WORLD,n, &x);</div><div>  ierr = VecCreateSeq(PETSC_COMM_WORLD,m, &residue);</div><div><br></div><div><br></div><div>  /* Set the solver type at the command-line */</div>
<div>  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);</div><div>  ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);</div><div>  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);</div><div>  ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);</div>
<div>  ierr = KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);</div><div>  /* Set runtime options, e.g.,</div><div>         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol></div>
<div>     These options will override those specified above as long as</div><div>     KSPSetFromOptions() is called _after_ any other customization</div><div>     routines */</div><div>  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);</div>
<div><br></div><div>  /* Initial guess for the krylov method set to zero*/</div><div>  PetscScalar p = 0;</div><div>  ierr = VecSet(x,p);CHKERRQ(ierr);</div><div>  ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);</div>
<div>  /* Solve linear system */</div><div>  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);</div><div>  /* View solver info; we could instead use the option -ksp_view to print this info to the screen at the conclusion of KSPSolve().*/</div>
<div>  ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);</div><div><br></div><div><br></div><div>  /*View the numerical solution and the residue vector*/</div><div>  printf("-------------------\nThe numerical solution is x = \n");</div>
<div>  ierr = VecView(x,PETSC_VIEWER_STDOUT_SELF);</div><div><br></div><div><br></div><div>  printf("-------------------\nThe residue vector Ax - b is = \n");</div><div>  ierr = MatMult(A,x,residue);CHKERRQ(ierr);</div>
<div>  ierr = VecAXPY(residue,-1.0,b);CHKERRQ(ierr); </div><div><br></div><div>  ierr = VecView(residue,PETSC_VIEWER_STDOUT_SELF);</div><div><br></div><div>  /* Clean up */</div><div>  ierr = VecDestroy(&x);CHKERRQ(ierr); </div>
<div>  ierr = VecDestroy(&residue);CHKERRQ(ierr);</div><div>  ierr = VecDestroy(&b);CHKERRQ(ierr); </div><div>  ierr = MatDestroy(&A);CHKERRQ(ierr);</div><div>  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);</div>
<div><br></div><div>  ierr = PetscFinalize();</div><div>  return 0;</div><div>}</div></div><div><br></div><div>  </div>