On Fri, Aug 3, 2012 at 5:38 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">







<div>
<font face="Calibri, sans-serif">
<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" target="_blank"><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" target="_blank"><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></font></div></blockquote><div><br></div><div>It is ILU, the default preconditioner, that requires a square matrix. Try -pc_type none</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"><div><font face="Calibri, sans-serif">
<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><span class="HOEnZb"><font color="#888888">
<div> </div>
<div>Shuangshuang</div>
<div><font face="Times New Roman, serif"> </font></div>
<div><font face="Times New Roman, serif"> </font></div>
</font></span></font>
</div>

</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>