[petsc-users] error when using ksplsqr for overdetermined system

Jin, Shuangshuang Shuangshuang.Jin at pnnl.gov
Fri Aug 3 17:38:15 CDT 2012


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


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120803/b8271fd4/attachment-0001.html>


More information about the petsc-users mailing list