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