[petsc-users] loading a matrix in PETSc : Problems

[Gaurish Telang]

Hi,


Given matrix A and column vector b. I am trying to do a Least squares
problem Ax=b on PETSc.

A has dimenstion 2683x1274

I have A given to me in a text file in the form of three columns (Row,
Column,NonZeroElement).

E.g. If a matrix of dimension 2x3  such as

0 5 3
5 6 0
 would be given in the text file as:

1 2 5
1 3 3
2 1 5
2 2 6

The third column just lists the non-zero elements and the corresponding
first two columsn the position,.

Now I am trying to input this textfile data into  PETSc matrix A. But my
matrix takes forever to load.

The same code on small matrices of dimension say 4x3 works perfectly.

I have pasted my PETSc code below. Any suggestions on improving it will be
helpful.

***Explanation of my code****:

For my code I start with two text files: N1.mat and b1.mat---->

N1.mat containing the matrix A in (Row,Column,Nonzeroelement)  format and
b1.mat containing the right hand side b.

Then I read in the three columns of N1.mat as a matrix B which  I use to
re-create the matrix A. The dimension of B is lenB x 3 where lenB is the
number of lines in the text file containing B(39051 in my case).
b1.mat gets read in directly as a the vector b.

%=========================================================================
%=========================================================================
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;

int tempi,tempj;

Mat A;

FILE *fp;

MPI_Comm comm;

PetscInt n=2683,m=1274,index,lenB=39051;

PetscScalar scalar,rhs[n],B[lenB][3];

PetscErrorCode ierr;

PetscInt i,j;

comm = MPI_COMM_SELF;

/* 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);

/* 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");

if (fp==NULL)

{

fprintf(stderr, "Cannot open file");

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);

}

}

index=0;

/*Putting x into final form */

for (i=0; i<n; ++i)

{

ierr= VecSetValues(b,1,&index,&rhs[i],INSERT_VALUES);CHKERRQ(ierr); /* One
insertion per step. Thats what the 1 in second argument stands for */

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);

/*
------------------------------------------------------------------------------------------------------------------------------------------------------
*/

ierr=MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);

ierr = MatSetSizes(A,n,m,n,m);CHKERRQ(ierr); /* Setting the matrix size */

ierr = MatSetFromOptions(A);CHKERRQ(ierr); /* Setting the matrix type
(shared, mpi etc) */

/* Reading in the matrix from the file matrix.txt */

fp=fopen("N1.mat","r");

if (fp==NULL)

{

fprintf(stderr, "Cannot open file");

exit(1);

}

for (i = 0; i < lenB; i++)

{

for (j = 0; j < 3; j++)

{

if (fscanf(fp,"%lf", &B[i][j]) != 1)

{

fprintf(stderr, "Failed to read matrix[%d][%d]\n", i, j);

exit(1);

}

}

}

/* Initializing the matrix A to zero matrix */

for(i=0;i<n;++i)

{

for(j=0;j<m;++j)

{

scalar=0.0;

ierr=MatSetValues(A,1,&i,1,&j,&scalar,INSERT_VALUES);CHKERRQ(ierr);

}

}

ierr=MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

ierr=MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

/* Inserting the non-zero elements. */

for(i=0;i<lenB;++i)

{

tempi=(int)B[i][0]-1;

tempj=(int)B[i][1]-1;

ierr=MatSetValues(A,1,&tempi,1,&tempj,&B[i][2],INSERT_VALUES);CHKERRQ(ierr);

}

ierr=MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

ierr=MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

//ierr=MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

/*

Destroy any objects created.

*/

ierr=VecDestroy(b);CHKERRQ(ierr);

ierr=MatDestroy(A);CHKERRQ(ierr);

ierr=PetscFinalize();CHKERRQ(ierr);

return 0;

}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110116/f0c62de9/attachment.htm>