<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Jan 23, 2013 at 12:28 AM, Gaurish Telang <span dir="ltr"><<a href="mailto:gaurish108@gmail.com" target="_blank">gaurish108@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">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>
</blockquote><div><br></div><div style>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.<br>
<br>src/ksp/ksp/examples/tutorials/ex10.c automatically loads the matrix in parallel and solves non-square systems (use -ksp_type lsqr).</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<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>
</blockquote></div><br></div></div>