<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Exchange Server">
<!-- converted from rtf -->
<style><!-- .EmailQuote { margin-left: 1pt; padding-left: 4pt; border-left: #800000 2px solid; } --></style>
</head>
<body>
<font face="Calibri, sans-serif" size="2">
<div>Hi, all, I’m trying to use ksplsqr to solve an overdetermined system. </div>
<div>I have read the following threads regarding the topic:</div>
<div><font face="Times New Roman, serif"><a href="http://lists.mcs.anl.gov/pipermail/petsc-users/2011-January/007763.html"><font face="Calibri, sans-serif" color="#0000FF"><u>http://lists.mcs.anl.gov/pipermail/petsc-users/2011-January/007763.html</u></font></a></font></div>
<div><font face="Times New Roman, serif"><a href="http://lists.mcs.anl.gov/pipermail/petsc-users/2010-August/006784.html"><font face="Calibri, sans-serif" color="#0000FF"><u>http://lists.mcs.anl.gov/pipermail/petsc-users/2010-August/006784.html</u></font></a></font></div>
<div><font face="Times New Roman, serif"> </font></div>
<div>I defined a 4*3 matrix A and a vector b in sparse CSR format:</div>
<div>irow: 1 3 5 6 9</div>
<div>icol: 1 2 1 3 2 1 2 3</div>
<div>ival: 1 2 2 5 3 1 4 1</div>
<div>rhs: 1 2 3 4</div>
<div><font face="Times New Roman, serif"> </font></div>
<div>A = 1 2 0</div>
<div> 2 0 5</div>
<div> 0 3 0</div>
<div> 1 4 1</div>
<div>b = 1</div>
<div> 2</div>
<div> 3</div>
<div> 4</div>
<div><font face="Times New Roman, serif"> </font></div>
<div>I wrote the following code to solve this system:</div>
<div><font face="Times New Roman, serif"> </font></div>
<div>static char help[]="Reading in a matrix\n";</div>
<div> </div>
<div>#include<stdio.h></div>
<div>#include<math.h></div>
<div>#include<stdlib.h></div>
<div>#include "petscvec.h" /* This enables us to use vectors. */</div>
<div>#include "petscmat.h" /* This enables us to use Matrices. It includes the petscvec header file*/</div>
<div>#include "petscksp.h" /* Now we can solve linear systems. Solvers used are KSP. */</div>
<div> </div>
<div>int main(int argc, char **argv)</div>
<div>{</div>
<div> /* Declaration of Matrix A and some vectors x*/</div>
<div> Vec b,x;</div>
<div> Mat A;</div>
<div> FILE *fp;</div>
<div> </div>
<div> MPI_Comm comm;</div>
<div> PetscInt n=4,m=3,nz=8,index;</div>
<div> PetscErrorCode ierr;</div>
<div> PetscInt i,j;</div>
<div> comm = MPI_COMM_SELF;</div>
<div> PetscInt irow[n+1],icol[nz]; </div>
<div> PetscScalar ival[nz],rhs[n];</div>
<div> PetscInt *cnt,*col;</div>
<div> </div>
<div> KSP ksp;</div>
<div> </div>
<div> /* This part is needed for the help flag supplied at run-time*/</div>
<div> ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);</div>
<div> ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);</div>
<div> </div>
<div>/*===================================================================================*/</div>
<div> </div>
<div> fp=fopen("irow.txt","r");</div>
<div> if (fp==NULL)</div>
<div> {</div>
<div> fprintf(stderr, "Cannot open file\n");</div>
<div> exit(1);</div>
<div> }</div>
<div> </div>
<div> for (i = 0; i < n+1; i++)</div>
<div> {</div>
<div> if (fscanf(fp,"%d", &irow[i]) != 1)</div>
<div> {</div>
<div> fprintf(stderr, "Failed to read irow vector[%d]\n", i);</div>
<div> exit(1);</div>
<div> } </div>
<div> }</div>
<div> fclose(fp);</div>
<div> for (i = 0; i < n+1; i++) printf("irow[%d]=%d\n",i,irow[i]);</div>
<div> </div>
<div> fp=fopen("icol.txt","r");</div>
<div> if (fp==NULL)</div>
<div> {</div>
<div> fprintf(stderr, "Cannot open file\n");</div>
<div> exit(1);</div>
<div> }</div>
<div> </div>
<div> for (i = 0; i < nz; i++)</div>
<div> {</div>
<div> if (fscanf(fp,"%d", &icol[i]) != 1)</div>
<div> {</div>
<div> fprintf(stderr, "Failed to read icol vector[%d]\n", i);</div>
<div> exit(1);</div>
<div> }</div>
<div> }</div>
<div> fclose(fp);</div>
<div> for (i = 0; i < nz; i++) printf("icol[%d]=%d\n",i,icol[i]);</div>
<div> </div>
<div> fp=fopen("ival.txt","r");</div>
<div> if (fp==NULL)</div>
<div> {</div>
<div> fprintf(stderr, "Cannot open file\n");</div>
<div> exit(1);</div>
<div> }</div>
<div> </div>
<div> for (i = 0; i < nz; i++)</div>
<div> {</div>
<div> if (fscanf(fp,"%lf", &ival[i]) != 1)</div>
<div> {</div>
<div> fprintf(stderr, "Failed to read ival vector[%d]\n", i);</div>
<div> exit(1);</div>
<div> }</div>
<div> }</div>
<div> fclose(fp);</div>
<div> for (i = 0; i < nz; i++) printf("ival[%d]=%lf\n",i,ival[i]);</div>
<div> </div>
<div>/*===================================================================================*/</div>
<div> </div>
<div> // determine number of nonzeros per row in the new matrix</div>
<div> PetscMalloc((n+1)*sizeof(PetscInt),&cnt);</div>
<div> for (i=0;i<n;i++) {</div>
<div> cnt[i]=irow[i+1]-irow[i];</div>
<div> }</div>
<div> MatCreateSeqAIJ(PETSC_COMM_SELF,n,m,0,cnt,&A);</div>
<div> MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);</div>
<div> // copy over the matrix entries from the inputs</div>
<div> for (i=0;i<n;i++) {</div>
<div> PetscMalloc(cnt[i]*sizeof(PetscInt),&col);</div>
<div> for (j=0;j<cnt[i];j++) {</div>
<div> col[j]=icol[irow[i]+j-1]-1;</div>
<div> }</div>
<div> for (j=irow[i]-1;j<irow[i]-1+cnt[i];j++) </div>
<div> MatSetValues(A,1,&i,cnt[i],col,&ival[irow[i]-1],INSERT_VALUES);</div>
<div> } </div>
<div> MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);</div>
<div> MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);</div>
<div> ierr=PetscPrintf(PETSC_COMM_WORLD,"Matrix A:\n");CHKERRQ(ierr);</div>
<div> MatView(A,PETSC_VIEWER_STDOUT_SELF);</div>
<div> </div>
<div>/*===================================================================================*/</div>
<div> </div>
<div> /* Use options from the terminal to create a vector that is type shared or mpi. */</div>
<div> ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr); /* Vector creation */</div>
<div> ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRQ(ierr); /* Setting the vector size */</div>
<div> ierr = VecSetFromOptions(b);CHKERRQ(ierr); /* Setting the vector type (shared, mpi etc) */</div>
<div> /* The second argument is a PETSc scalar value.*/</div>
<div> ierr = VecSet(b,0);CHKERRQ(ierr);</div>
<div> </div>
<div> /* Reading in the RHS vector. */</div>
<div> /* Reading in the matrix from the file matrix.txt */</div>
<div> //fp=fopen("b1.mat","r");</div>
<div> fp=fopen("rhs.txt","r");</div>
<div> if (fp==NULL)</div>
<div> {</div>
<div> fprintf(stderr, "Cannot open file\n");</div>
<div> exit(1);</div>
<div> }</div>
<div> </div>
<div> for (i = 0; i < n; i++)</div>
<div> {</div>
<div> if (fscanf(fp,"%lf", &rhs[i]) != 1)</div>
<div> {</div>
<div> fprintf(stderr, "Failed to read rhs vector[%d]\n", i);</div>
<div> exit(1);</div>
<div> }</div>
<div> }</div>
<div> fclose(fp);</div>
<div> </div>
<div> index=0;</div>
<div> /*Putting x into final form */</div>
<div> for (i=0; i<n; ++i)</div>
<div> {</div>
<div> /* Oneinsertion per step. Thats what the 1 in second argument stands for */</div>
<div> ierr= VecSetValues(b,1,&index,&rhs[i],INSERT_VALUES);CHKERRQ(ierr); </div>
<div> index=index+1;</div>
<div> } /* The third and fourth arguments are addresses. The fifth argument is IORA */</div>
<div> </div>
<div> /* Assembling the vector. */</div>
<div> ierr= VecAssemblyBegin(b);CHKERRQ(ierr);</div>
<div> ierr=VecAssemblyEnd(b);CHKERRQ(ierr);</div>
<div> </div>
<div> /* Viewing the changed vector. */</div>
<div> ierr=PetscPrintf(PETSC_COMM_WORLD,"Vector b:\n");CHKERRQ(ierr);</div>
<div> ierr=VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);</div>
<div> </div>
<div>/*===================================================================================*/</div>
<div> </div>
<div> /* solve kspslqr system */</div>
<div> VecCreate(PETSC_COMM_WORLD,&x);</div>
<div> VecSetSizes(x,PETSC_DECIDE,m);</div>
<div> VecSetFromOptions(x);</div>
<div> </div>
<div> KSPCreate(PETSC_COMM_WORLD, &ksp); </div>
<div> KSPSetType(ksp, KSPLSQR);</div>
<div> KSPSetTolerances(ksp,1.0e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);</div>
<div> KSPSetOperators(ksp, A, A, SAME_PRECONDITIONER);</div>
<div> KSPSetFromOptions(ksp);</div>
<div> KSPSolve(ksp, b, x);</div>
<div> </div>
<div> PetscInt num_iters;</div>
<div> KSPGetIterationNumber(ksp, &num_iters); CHKERRQ(ierr);</div>
<div> </div>
<div> PetscReal rnorm;</div>
<div> ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr);</div>
<div> </div>
<div> KSPConvergedReason reason;</div>
<div> ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr);</div>
<div> </div>
<div> printf("KSPGetIterationNumber %d\n",num_iters);</div>
<div> printf("KSPGetResidualNorm %f\n",rnorm);</div>
<div> printf("KSPConvergedReason %d\n",reason);</div>
<div> </div>
<div> //VecView(x,PETSC_VIEWER_STDOUT_WORLD);</div>
<div> //PetscViewerASCIIOpen(PETSC_COMM_WORLD, "x.data", &viewer);</div>
<div> VecView(x,PETSC_VIEWER_STDOUT_WORLD);</div>
<div> KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);</div>
<div> </div>
<div> /*</div>
<div> Destroy any objects created.</div>
<div> */</div>
<div> ierr=VecDestroy(&b);CHKERRQ(ierr);</div>
<div> ierr=MatDestroy(&A);CHKERRQ(ierr);</div>
<div> ierr=PetscFinalize();CHKERRQ(ierr);</div>
<div> </div>
<div> return 0;</div>
<div>}</div>
<div><font face="Times New Roman, serif"> </font></div>
<div>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:</div>
<div><font face="Times New Roman, serif"> </font></div>
<div>[d3m956@olympus tutorials]$ ./ss</div>
<div>irow[0]=1</div>
<div>irow[1]=3</div>
<div>irow[2]=5</div>
<div>irow[3]=6</div>
<div>irow[4]=9</div>
<div>icol[0]=1</div>
<div>icol[1]=2</div>
<div>icol[2]=1</div>
<div>icol[3]=3</div>
<div>icol[4]=2</div>
<div>icol[5]=1</div>
<div>icol[6]=2</div>
<div>icol[7]=3</div>
<div>ival[0]=1.000000</div>
<div>ival[1]=2.000000</div>
<div>ival[2]=2.000000</div>
<div>ival[3]=5.000000</div>
<div>ival[4]=3.000000</div>
<div>ival[5]=1.000000</div>
<div>ival[6]=4.000000</div>
<div>ival[7]=1.000000</div>
<div>Matrix A:</div>
<div>Matrix Object: 1 MPI processes</div>
<div> type: seqaij</div>
<div>row 0: (0, 1) (1, 2)</div>
<div>row 1: (0, 2) (2, 5)</div>
<div>row 2: (1, 3)</div>
<div>row 3: (0, 1) (1, 4) (2, 1)</div>
<div>Vector b:</div>
<div>Vector Object: 1 MPI processes</div>
<div> type: seq</div>
<div>1</div>
<div>2</div>
<div>3</div>
<div>4</div>
<div>[0]PETSC ERROR: --------------------- Error Message ------------------------------------</div>
<div>[0]PETSC ERROR: Invalid argument!</div>
<div>[0]PETSC ERROR: Must be square matrix, rows 4 columns 3!</div>
<div>[0]PETSC ERROR: ------------------------------------------------------------------------</div>
<div>[0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 0, Tue Jun 5 14:20:42 CDT 2012</div>
<div>[0]PETSC ERROR: See docs/changes/index.html for recent updates.</div>
<div>[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.</div>
<div>[0]PETSC ERROR: See docs/index.html for manual pages.</div>
<div>[0]PETSC ERROR: ------------------------------------------------------------------------</div>
<div>[0]PETSC ERROR: ./ss on a arch-linu named olympus.local by d3m956 Fri Aug 3 15:34:07 2012</div>
<div>[0]PETSC ERROR: Libraries linked from /pic/projects/mca/ss/PETSC/petsc-3.3-p0/arch-linux2-c-debug/lib</div>
<div>[0]PETSC ERROR: Configure run at Thu Jun 14 17:00:19 2012</div>
<div>[0]PETSC ERROR: Configure options --with-cc=gcc --with-fc=gfortran --download-f-blas-lapack --download-mpich</div>
<div>[0]PETSC ERROR: ------------------------------------------------------------------------</div>
<div>[0]PETSC ERROR: MatGetOrdering() line 228 in src/mat/order/sorder.c</div>
<div>[0]PETSC ERROR: PCSetUp_ILU() line 206 in src/ksp/pc/impls/factor/ilu/ilu.c</div>
<div>[0]PETSC ERROR: PCSetUp() line 832 in src/ksp/pc/interface/precon.c</div>
<div>[0]PETSC ERROR: KSPSetUp() line 278 in src/ksp/ksp/interface/itfunc.c</div>
<div>[0]PETSC ERROR: KSPSolve() line 402 in src/ksp/ksp/interface/itfunc.c</div>
<div>KSPGetIterationNumber 0</div>
<div>KSPGetResidualNorm 0.000000</div>
<div>KSPConvergedReason 0</div>
<div>Vector Object: 1 MPI processes</div>
<div> type: seq</div>
<div>0</div>
<div>0</div>
<div>0</div>
<div>KSP Object: 1 MPI processes</div>
<div> type: lsqr</div>
<div> maximum iterations=10000, initial guess is zero</div>
<div> tolerances: relative=1e-06, absolute=1e-50, divergence=10000</div>
<div> left preconditioning</div>
<div> using UNPRECONDITIONED norm type for convergence test</div>
<div>PC Object: 1 MPI processes</div>
<div> type: ilu</div>
<div> ILU: out-of-place factorization</div>
<div> 0 levels of fill</div>
<div> tolerance for zero pivot 2.22045e-14</div>
<div> using diagonal shift to prevent zero pivot</div>
<div> matrix ordering: natural</div>
<div> linear system matrix = precond matrix:</div>
<div> Matrix Object: 1 MPI processes</div>
<div> type: seqaij</div>
<div> rows=4, cols=3</div>
<div> total: nonzeros=8, allocated nonzeros=8</div>
<div> total number of mallocs used during MatSetValues calls =0</div>
<div> not using I-node routines</div>
<div><font face="Times New Roman, serif"> </font></div>
<div>I don’t understand why the lsqr solver requires a square matrix “Must be square matrix, rows 4 columns 3!”</div>
<div><font face="Times New Roman, serif"> </font></div>
<div>Can anyone look into the problem for me?</div>
<div> </div>
<div>Thanks a lot!</div>
<div> </div>
<div>Shuangshuang</div>
<div><font face="Times New Roman, serif"> </font></div>
<div><font face="Times New Roman, serif"> </font></div>
</font>
</body>
</html>