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/20120803/7940e11d/attachment.html>
