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

[Matthew Knepley]

You need to preallocate the matrix memory. I would read through the file
once, count the number
of nonzeros in each row, preallocate the matrix, then read the file again to
put in the data.

   Matt

On Sun, Jan 16, 2011 at 3:21 PM, Gaurish Telang <gaurish108 at gmail.com>wrote:

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



-- 
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/20110116/1890a489/attachment.htm>