On Sun, Aug 5, 2012 at 3:49 PM, Jin, Shuangshuang <span dir="ltr"><<a href="mailto:Shuangshuang.Jin@pnnl.gov" target="_blank">Shuangshuang.Jin@pnnl.gov</a>></span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi, Mat, I replaced the inputs with a larger system (the A matrix has 27090 rows and 14959 columns). With the -pc_type none set up as the runtime option, the system can be solved correctly, but the new problem is the solving time is very slow, taking about 75 seconds to run KSPSolve for the lsqr problem. Is there any faster preconditioning type we can use here instead of "none"?<br>
</blockquote><div><br></div><div>I have no idea what system you are solving. I suggest consulting the literature first. If you know what</div><div>preconditioner you want, we can try to help.</div><div><br></div><div>About parallelism, <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html#computers">http://www.mcs.anl.gov/petsc/documentation/faq.html#computers</a></div>
<div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
[d3m956@olympus tutorials]$ mpiexec -n 1 ss -pc_type none<br>
Pre KSPSolve Elapsed Time: 0.195198<br>
KSPSolve Elapsed Time: 74.788414<br>
KSPGetIterationNumber 10000<br>
KSPGetResidualNorm 17921.356837<br>
KSPConvergedReason -3<br>
KSP Object: 1 MPI processes<br>
type: lsqr<br>
maximum iterations=10000, initial guess is zero<br>
tolerances: relative=1e-06, absolute=1e-50, divergence=10000<br>
left preconditioning<br>
using UNPRECONDITIONED norm type for convergence test<br>
PC Object: 1 MPI processes<br>
type: none<br>
linear system matrix = precond matrix:<br>
Matrix Object: 1 MPI processes<br>
type: seqaij<br>
rows=27090, cols=14959<br>
total: nonzeros=70903, allocated nonzeros=70903<br>
total number of mallocs used during MatSetValues calls =0<br>
not using I-node routines<br>
Overall Elapsed Time: 75.155196<br>
<br>
Another question is how to run lsqr in parallel. I tried to use 4 processors for example, but it didn't help to speed up.<br>
<br>
[d3m956@olympus tutorials]$ mpiexec -n 4 ss -pc_type none<br>
Pre KSPSolve Elapsed Time: 0.177819<br>
Pre KSPSolve Elapsed Time: 0.180807<br>
Pre KSPSolve Elapsed Time: 0.179783<br>
Pre KSPSolve Elapsed Time: 0.181975<br>
KSPSolve Elapsed Time: 73.107033<br>
KSPGetIterationNumber 10000<br>
KSPGetResidualNorm 17921.356837<br>
KSPConvergedReason -3<br>
KSP Object: 1 MPI processes<br>
type: lsqr<br>
maximum iterations=10000, initial guess is zero<br>
tolerances: relative=1e-06, absolute=1e-50, divergence=10000<br>
left preconditioning<br>
using UNPRECONDITIONED norm type for convergence test<br>
PC Object: 1 MPI processes<br>
type: none<br>
linear system matrix = precond matrix:<br>
Matrix Object: 1 MPI processes<br>
type: seqaij<br>
rows=27090, cols=14959<br>
total: nonzeros=70903, allocated nonzeros=70903<br>
total number of mallocs used during MatSetValues calls =0<br>
not using I-node routines<br>
Overall Elapsed Time: 73.455209<br>
KSPSolve Elapsed Time: 73.425921<br>
KSPGetIterationNumber 10000<br>
KSPGetResidualNorm 17921.356837<br>
KSPConvergedReason -3<br>
KSP Object: 1 MPI processes<br>
type: lsqr<br>
maximum iterations=10000, initial guess is zero<br>
tolerances: relative=1e-06, absolute=1e-50, divergence=10000<br>
left preconditioning<br>
using UNPRECONDITIONED norm type for convergence test<br>
PC Object: 1 MPI processes<br>
type: none<br>
linear system matrix = precond matrix:<br>
Matrix Object: 1 MPI processes<br>
type: seqaij<br>
rows=27090, cols=14959<br>
total: nonzeros=70903, allocated nonzeros=70903<br>
total number of mallocs used during MatSetValues calls =0<br>
not using I-node routines<br>
Overall Elapsed Time: 73.771948<br>
KSPSolve Elapsed Time: 76.629414<br>
KSPGetIterationNumber 10000<br>
KSPGetResidualNorm 17921.356837<br>
KSPConvergedReason -3<br>
KSPSolve Elapsed Time: 76.698510<br>
KSPGetIterationNumber 10000<br>
KSPGetResidualNorm 17921.356837<br>
KSPConvergedReason -3<br>
KSP Object: 1 MPI processes<br>
type: lsqr<br>
maximum iterations=10000, initial guess is zero<br>
tolerances: relative=1e-06, absolute=1e-50, divergence=10000<br>
left preconditioning<br>
using UNPRECONDITIONED norm type for convergence test<br>
PC Object: 1 MPI processes<br>
type: none<br>
linear system matrix = precond matrix:<br>
Matrix Object: 1 MPI processes<br>
type: seqaij<br>
rows=27090, cols=14959<br>
total: nonzeros=70903, allocated nonzeros=70903<br>
total number of mallocs used during MatSetValues calls =0<br>
not using I-node routines<br>
Overall Elapsed Time: 77.174527<br>
KSP Object: 1 MPI processes<br>
type: lsqr<br>
maximum iterations=10000, initial guess is zero<br>
tolerances: relative=1e-06, absolute=1e-50, divergence=10000<br>
left preconditioning<br>
using UNPRECONDITIONED norm type for convergence test<br>
PC Object: 1 MPI processes<br>
type: none<br>
linear system matrix = precond matrix:<br>
Matrix Object: 1 MPI processes<br>
type: seqaij<br>
rows=27090, cols=14959<br>
total: nonzeros=70903, allocated nonzeros=70903<br>
total number of mallocs used during MatSetValues calls =0<br>
not using I-node routines<br>
Overall Elapsed Time: 77.247951<br>
<br>
Is there anything I need to do in the code to parallelize the ksp lsqr solver other than specify the mpiexec -n <NumOfPrcessors> option in the command line?<br>
<br>
Thanks,<br>
Shuangshuang<br>
<br>
<br>
On 8/3/12 4:35 PM, "Matthew Knepley" <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
<br>
On Fri, Aug 3, 2012 at 5:38 PM, Jin, Shuangshuang <<a href="mailto:Shuangshuang.Jin@pnnl.gov">Shuangshuang.Jin@pnnl.gov</a>> wrote:<br>
Hi, all, I'm trying to use ksplsqr to solve an overdetermined system.<br>
I have read the following threads regarding the topic:<br>
<a href="http://lists.mcs.anl.gov/pipermail/petsc-users/2011-January/007763.html" target="_blank">http://lists.mcs.anl.gov/pipermail/petsc-users/2011-January/007763.html</a> <<a href="http://lists.mcs.anl.gov/pipermail/petsc-users/2011-January/007763.html" target="_blank">http://lists.mcs.anl.gov/pipermail/petsc-users/2011-January/007763.html</a>><br>
<a href="http://lists.mcs.anl.gov/pipermail/petsc-users/2010-August/006784.html" target="_blank">http://lists.mcs.anl.gov/pipermail/petsc-users/2010-August/006784.html</a> <<a href="http://lists.mcs.anl.gov/pipermail/petsc-users/2010-August/006784.html" target="_blank">http://lists.mcs.anl.gov/pipermail/petsc-users/2010-August/006784.html</a>><br>
<br>
I defined a 4*3 matrix A and a vector b in sparse CSR format:<br>
irow: 1 3 5 6 9<br>
icol: 1 2 1 3 2 1 2 3<br>
ival: 1 2 2 5 3 1 4 1<br>
rhs: 1 2 3 4<br>
<br>
A = 1 2 0<br>
2 0 5<br>
0 3 0<br>
1 4 1<br>
b = 1<br>
2<br>
3<br>
4<br>
<br>
I wrote the following code to solve this system:<br>
<br>
static char help[]="Reading in a matrix\n";<br>
<br>
#include<stdio.h><br>
#include<math.h><br>
#include<stdlib.h><br>
#include "petscvec.h" /* This enables us to use vectors. */<br>
#include "petscmat.h" /* This enables us to use Matrices. It includes the petscvec header file*/<br>
#include "petscksp.h" /* Now we can solve linear systems. Solvers used are KSP. */<br>
<br>
int main(int argc, char **argv)<br>
{<br>
/* Declaration of Matrix A and some vectors x*/<br>
Vec b,x;<br>
Mat A;<br>
FILE *fp;<br>
<br>
MPI_Comm comm;<br>
PetscInt n=4,m=3,nz=8,index;<br>
PetscErrorCode ierr;<br>
PetscInt i,j;<br>
comm = MPI_COMM_SELF;<br>
PetscInt irow[n+1],icol[nz];<br>
PetscScalar ival[nz],rhs[n];<br>
PetscInt *cnt,*col;<br>
<br>
KSP ksp;<br>
<br>
/* This part is needed for the help flag supplied at run-time*/<br>
ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);<br>
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);<br>
<br>
/*===================================================================================*/<br>
<br>
fp=fopen("irow.txt","r");<br>
if (fp==NULL)<br>
{<br>
fprintf(stderr, "Cannot open file\n");<br>
exit(1);<br>
}<br>
<br>
for (i = 0; i < n+1; i++)<br>
{<br>
if (fscanf(fp,"%d", &irow[i]) != 1)<br>
{<br>
fprintf(stderr, "Failed to read irow vector[%d]\n", i);<br>
exit(1);<br>
}<br>
}<br>
fclose(fp);<br>
for (i = 0; i < n+1; i++) printf("irow[%d]=%d\n",i,irow[i]);<br>
<br>
fp=fopen("icol.txt","r");<br>
if (fp==NULL)<br>
{<br>
fprintf(stderr, "Cannot open file\n");<br>
exit(1);<br>
}<br>
<br>
for (i = 0; i < nz; i++)<br>
{<br>
if (fscanf(fp,"%d", &icol[i]) != 1)<br>
{<br>
fprintf(stderr, "Failed to read icol vector[%d]\n", i);<br>
exit(1);<br>
}<br>
}<br>
fclose(fp);<br>
for (i = 0; i < nz; i++) printf("icol[%d]=%d\n",i,icol[i]);<br>
<br>
fp=fopen("ival.txt","r");<br>
if (fp==NULL)<br>
{<br>
fprintf(stderr, "Cannot open file\n");<br>
exit(1);<br>
}<br>
<br>
for (i = 0; i < nz; i++)<br>
{<br>
if (fscanf(fp,"%lf", &ival[i]) != 1)<br>
{<br>
fprintf(stderr, "Failed to read ival vector[%d]\n", i);<br>
exit(1);<br>
}<br>
}<br>
fclose(fp);<br>
for (i = 0; i < nz; i++) printf("ival[%d]=%lf\n",i,ival[i]);<br>
<br>
/*===================================================================================*/<br>
<br>
// determine number of nonzeros per row in the new matrix<br>
PetscMalloc((n+1)*sizeof(PetscInt),&cnt);<br>
for (i=0;i<n;i++) {<br>
cnt[i]=irow[i+1]-irow[i];<br>
}<br>
MatCreateSeqAIJ(PETSC_COMM_SELF,n,m,0,cnt,&A);<br>
MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);<br>
// copy over the matrix entries from the inputs<br>
for (i=0;i<n;i++) {<br>
PetscMalloc(cnt[i]*sizeof(PetscInt),&col);<br>
for (j=0;j<cnt[i];j++) {<br>
col[j]=icol[irow[i]+j-1]-1;<br>
}<br>
for (j=irow[i]-1;j<irow[i]-1+cnt[i];j++)<br>
MatSetValues(A,1,&i,cnt[i],col,&ival[irow[i]-1],INSERT_VALUES);<br>
}<br>
MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);<br>
MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);<br>
ierr=PetscPrintf(PETSC_COMM_WORLD,"Matrix A:\n");CHKERRQ(ierr);<br>
MatView(A,PETSC_VIEWER_STDOUT_SELF);<br>
<br>
/*===================================================================================*/<br>
<br>
/* Use options from the terminal to create a vector that is type shared or mpi. */<br>
ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr); /* Vector creation */<br>
ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRQ(ierr); /* Setting the vector size */<br>
ierr = VecSetFromOptions(b);CHKERRQ(ierr); /* Setting the vector type (shared, mpi etc) */<br>
/* The second argument is a PETSc scalar value.*/<br>
ierr = VecSet(b,0);CHKERRQ(ierr);<br>
<br>
/* Reading in the RHS vector. */<br>
/* Reading in the matrix from the file matrix.txt */<br>
//fp=fopen("b1.mat","r");<br>
fp=fopen("rhs.txt","r");<br>
if (fp==NULL)<br>
{<br>
fprintf(stderr, "Cannot open file\n");<br>
exit(1);<br>
}<br>
<br>
for (i = 0; i < n; i++)<br>
{<br>
if (fscanf(fp,"%lf", &rhs[i]) != 1)<br>
{<br>
fprintf(stderr, "Failed to read rhs vector[%d]\n", i);<br>
exit(1);<br>
}<br>
}<br>
fclose(fp);<br>
<br>
index=0;<br>
/*Putting x into final form */<br>
for (i=0; i<n; ++i)<br>
{<br>
/* Oneinsertion per step. Thats what the 1 in second argument stands for */<br>
ierr= VecSetValues(b,1,&index,&rhs[i],INSERT_VALUES);CHKERRQ(ierr);<br>
index=index+1;<br>
} /* The third and fourth arguments are addresses. The fifth argument is IORA */<br>
<br>
/* Assembling the vector. */<br>
ierr= VecAssemblyBegin(b);CHKERRQ(ierr);<br>
ierr=VecAssemblyEnd(b);CHKERRQ(ierr);<br>
<br>
/* Viewing the changed vector. */<br>
ierr=PetscPrintf(PETSC_COMM_WORLD,"Vector b:\n");CHKERRQ(ierr);<br>
ierr=VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);<br>
<br>
/*===================================================================================*/<br>
<br>
/* solve kspslqr system */<br>
VecCreate(PETSC_COMM_WORLD,&x);<br>
VecSetSizes(x,PETSC_DECIDE,m);<br>
VecSetFromOptions(x);<br>
<br>
KSPCreate(PETSC_COMM_WORLD, &ksp);<br>
KSPSetType(ksp, KSPLSQR);<br>
KSPSetTolerances(ksp,1.0e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);<br>
KSPSetOperators(ksp, A, A, SAME_PRECONDITIONER);<br>
KSPSetFromOptions(ksp);<br>
KSPSolve(ksp, b, x);<br>
<br>
PetscInt num_iters;<br>
KSPGetIterationNumber(ksp, &num_iters); CHKERRQ(ierr);<br>
<br>
PetscReal rnorm;<br>
ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr);<br>
<br>
KSPConvergedReason reason;<br>
ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr);<br>
<br>
printf("KSPGetIterationNumber %d\n",num_iters);<br>
printf("KSPGetResidualNorm %f\n",rnorm);<br>
printf("KSPConvergedReason %d\n",reason);<br>
<br>
//VecView(x,PETSC_VIEWER_STDOUT_WORLD);<br>
//PetscViewerASCIIOpen(PETSC_COMM_WORLD, "x.data", &viewer);<br>
VecView(x,PETSC_VIEWER_STDOUT_WORLD);<br>
KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);<br>
<br>
/*<br>
Destroy any objects created.<br>
*/<br>
ierr=VecDestroy(&b);CHKERRQ(ierr);<br>
ierr=MatDestroy(&A);CHKERRQ(ierr);<br>
ierr=PetscFinalize();CHKERRQ(ierr);<br>
<br>
return 0;<br>
}<br>
<br>
This code read the inputs data and formed the A and b matrix correctly as I can view them in the output. However, I got "Invalid argument!" error as shown below:<br>
<br>
[d3m956@olympus tutorials]$ ./ss<br>
irow[0]=1<br>
irow[1]=3<br>
irow[2]=5<br>
irow[3]=6<br>
irow[4]=9<br>
icol[0]=1<br>
icol[1]=2<br>
icol[2]=1<br>
icol[3]=3<br>
icol[4]=2<br>
icol[5]=1<br>
icol[6]=2<br>
icol[7]=3<br>
ival[0]=1.000000<br>
ival[1]=2.000000<br>
ival[2]=2.000000<br>
ival[3]=5.000000<br>
ival[4]=3.000000<br>
ival[5]=1.000000<br>
ival[6]=4.000000<br>
ival[7]=1.000000<br>
Matrix A:<br>
Matrix Object: 1 MPI processes<br>
type: seqaij<br>
row 0: (0, 1) (1, 2)<br>
row 1: (0, 2) (2, 5)<br>
row 2: (1, 3)<br>
row 3: (0, 1) (1, 4) (2, 1)<br>
Vector b:<br>
Vector Object: 1 MPI processes<br>
type: seq<br>
1<br>
2<br>
3<br>
4<br>
[0]PETSC ERROR: --------------------- Error Message ------------------------------------<br>
[0]PETSC ERROR: Invalid argument!<br>
[0]PETSC ERROR: Must be square matrix, rows 4 columns 3!<br>
[0]PETSC ERROR: ------------------------------------------------------------------------<br>
[0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 0, Tue Jun 5 14:20:42 CDT 2012<br>
[0]PETSC ERROR: See docs/changes/index.html for recent updates.<br>
[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.<br>
[0]PETSC ERROR: See docs/index.html for manual pages.<br>
[0]PETSC ERROR: ------------------------------------------------------------------------<br>
[0]PETSC ERROR: ./ss on a arch-linu named olympus.local by d3m956 Fri Aug 3 15:34:07 2012<br>
[0]PETSC ERROR: Libraries linked from /pic/projects/mca/ss/PETSC/petsc-3.3-p0/arch-linux2-c-debug/lib<br>
[0]PETSC ERROR: Configure run at Thu Jun 14 17:00:19 2012<br>
[0]PETSC ERROR: Configure options --with-cc=gcc --with-fc=gfortran --download-f-blas-lapack --download-mpich<br>
[0]PETSC ERROR: ------------------------------------------------------------------------<br>
[0]PETSC ERROR: MatGetOrdering() line 228 in src/mat/order/sorder.c<br>
[0]PETSC ERROR: PCSetUp_ILU() line 206 in src/ksp/pc/impls/factor/ilu/ilu.c<br>
<br>
It is ILU, the default preconditioner, that requires a square matrix. Try -pc_type none<br>
<br>
Matt<br>
<br>
[0]PETSC ERROR: PCSetUp() line 832 in src/ksp/pc/interface/precon.c<br>
[0]PETSC ERROR: KSPSetUp() line 278 in src/ksp/ksp/interface/itfunc.c<br>
[0]PETSC ERROR: KSPSolve() line 402 in src/ksp/ksp/interface/itfunc.c<br>
KSPGetIterationNumber 0<br>
KSPGetResidualNorm 0.000000<br>
KSPConvergedReason 0<br>
Vector Object: 1 MPI processes<br>
type: seq<br>
0<br>
0<br>
0<br>
KSP Object: 1 MPI processes<br>
type: lsqr<br>
maximum iterations=10000, initial guess is zero<br>
tolerances: relative=1e-06, absolute=1e-50, divergence=10000<br>
left preconditioning<br>
using UNPRECONDITIONED norm type for convergence test<br>
PC Object: 1 MPI processes<br>
type: ilu<br>
ILU: out-of-place factorization<br>
0 levels of fill<br>
tolerance for zero pivot 2.22045e-14<br>
using diagonal shift to prevent zero pivot<br>
matrix ordering: natural<br>
linear system matrix = precond matrix:<br>
Matrix Object: 1 MPI processes<br>
type: seqaij<br>
rows=4, cols=3<br>
total: nonzeros=8, allocated nonzeros=8<br>
total number of mallocs used during MatSetValues calls =0<br>
not using I-node routines<br>
<br>
I don't understand why the lsqr solver requires a square matrix "Must be square matrix, rows 4 columns 3!"<br>
<br>
Can anyone look into the problem for me?<br>
<br>
Thanks a lot!<br>
<span class="HOEnZb"><font color="#888888"><br>
Shuangshuang<br>
<br>
<br>
<br>
<br>
</font></span></blockquote></div><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener<br>