You need to preallocate the matrix memory. I would read through the file once, count the number<div>of nonzeros in each row, preallocate the matrix, then read the file again to put in the data.</div><div><br></div><div> Matt<br>
<br><div class="gmail_quote">On Sun, Jan 16, 2011 at 3:21 PM, Gaurish Telang <span dir="ltr"><<a href="mailto:gaurish108@gmail.com">gaurish108@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div>Hi,</div>
<div> </div>
<div> </div>
<div>Given matrix A and column vector b. I am trying to do a Least squares problem Ax=b on PETSc. </div>
<div> </div>
<div>A has dimenstion 2683x1274</div>
<div> </div>
<div>I have A given to me in a text file in the form of three columns (Row, Column,NonZeroElement).</div>
<div> </div>
<div>E.g. If a matrix of dimension 2x3 such as </div>
<div> </div>
<div>0 5 3</div>
<div>5 6 0 </div>
<div> would be given in the text file as:</div>
<div> </div>
<div>1 2 5</div>
<div>1 3 3</div>
<div>2 1 5</div>
<div>2 2 6</div>
<div> </div>
<div>The third column just lists the non-zero elements and the corresponding first two columsn the position,.</div>
<div> </div>
<div>Now I am trying to input this textfile data into PETSc matrix A. But my matrix takes forever to load. </div>
<div> </div>
<div>The same code on small matrices of dimension say 4x3 works perfectly. </div>
<div> </div>
<div>I have pasted my PETSc code below. Any suggestions on improving it will be helpful.</div>
<div> </div>
<div>***Explanation of my code****: </div>
<div> </div>
<div>For my code I start with two text files: N1.mat and b1.mat----> </div>
<div> </div>
<div>N1.mat containing the matrix A in (Row,Column,Nonzeroelement) format and b1.mat containing the right hand side b. </div>
<div> </div>
<div>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). <br>
</div>
<div>b1.mat gets read in directly as a the vector b. </div>
<div> </div>
<div>%=========================================================================</div>
<div>%========================================================================= </div>
<div><span lang="EN">static char help[]="Reading in a matrix\n";
<p>#include<stdio.h></p>
<p>#include<math.h></p>
<p>#include<stdlib.h></p>
<p>#include "petscvec.h" /* This enables us to use vectors. */</p>
<p>#include "petscmat.h" /* This enables us to use Matrices. It includes the petscvec header file*/</p>
<p>#include "petscksp.h" /* Now we can solve linear systems. Solvers used are KSP. */</p>
<p>int main(int argc, char **argv)</p>
<p>{</p>
<p>/* Declaration of Matrix A and some vectors x*/</p>
<p>Vec b; </p>
<p>int tempi,tempj;</p>
<p>Mat A;</p>
<p>FILE *fp;</p>
<p>MPI_Comm comm;</p>
<p></p>
<p>PetscInt n=2683,m=1274,index,lenB=39051;</p>
<p>PetscScalar scalar,rhs[n],B[lenB][3]; </p>
<p>PetscErrorCode ierr;</p>
<p></p>
<p>PetscInt i,j;</p>
<p>comm = MPI_COMM_SELF;</p>
<p>/* This part is needed for the help flag supplied at run-time*/</p>
<p>ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr); </p>
<p>ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);</p>
<p>/* Use options from the terminal to create a vector that is type shared or mpi. */</p>
<p>ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr); /* Vector creation */</p>
<p>ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRQ(ierr); /* Setting the vector size */</p>
<p>ierr = VecSetFromOptions(b);CHKERRQ(ierr); /* Setting the vector type (shared, mpi etc) */</p>
<p>/* The second argument is a PETSc scalar value.*/</p>
<p>ierr = VecSet(b,0);CHKERRQ(ierr);</p>
<p>/* Reading in the RHS vector. */</p>
<p>/* Reading in the matrix from the file matrix.txt */</p>
<p>fp=fopen("b1.mat","r");</p>
<p>if (fp==NULL)</p>
<p>{</p>
<p>fprintf(stderr, "Cannot open file");</p>
<p>exit(1);</p>
<p>}</p>
<p>for (i = 0; i < n; i++)</p>
<p>{</p>
<p></p>
<p>if (fscanf(fp,"%lf", &rhs[i]) != 1)</p>
<p>{</p>
<p>fprintf(stderr, "Failed to read rhs vector[%d]\n", i);</p>
<p>exit(1);</p>
<p>}</p>
<p>}</p>
<p>index=0;</p>
<p>/*Putting x into final form */</p>
<p>for (i=0; i<n; ++i) </p>
<p>{</p>
<p>ierr= VecSetValues(b,1,&index,&rhs[i],INSERT_VALUES);CHKERRQ(ierr); /* One insertion per step. Thats what the 1 in second argument stands for */</p>
<p>index=index+1;</p>
<p>} /* The third and fourth arguments are addresses. The fifth argument is IORA */</p>
<p>/* Assembling the vector. */</p>
<p>ierr= VecAssemblyBegin(b);CHKERRQ(ierr);</p>
<p>ierr=VecAssemblyEnd(b);CHKERRQ(ierr);</p>
<p>/* Viewing the changed vector. */</p>
<p>ierr=PetscPrintf(PETSC_COMM_WORLD,"Vector b:\n");CHKERRQ(ierr);</p>
<p>ierr=VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); </p>
<p>/* ------------------------------------------------------------------------------------------------------------------------------------------------------ */</p>
<p>ierr=MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);</p>
<p>ierr = MatSetSizes(A,n,m,n,m);CHKERRQ(ierr); /* Setting the matrix size */</p>
<p>ierr = MatSetFromOptions(A);CHKERRQ(ierr); /* Setting the matrix type (shared, mpi etc) */</p>
<p>/* Reading in the matrix from the file matrix.txt */</p>
<p>fp=fopen("N1.mat","r");</p>
<p>if (fp==NULL)</p>
<p>{</p>
<p>fprintf(stderr, "Cannot open file");</p>
<p>exit(1);</p>
<p>}</p>
<p>for (i = 0; i < lenB; i++)</p>
<p>{</p>
<p>for (j = 0; j < 3; j++)</p>
<p>{</p>
<p>if (fscanf(fp,"%lf", &B[i][j]) != 1)</p>
<p>{</p>
<p>fprintf(stderr, "Failed to read matrix[%d][%d]\n", i, j);</p>
<p>exit(1);</p>
<p>}</p>
<p>}</p>
<p>}</p>
<p>/* Initializing the matrix A to zero matrix */</p>
<p>for(i=0;i<n;++i)</p>
<p>{</p>
<p>for(j=0;j<m;++j)</p>
<p>{</p>
<p>scalar=0.0;</p>
<p>ierr=MatSetValues(A,1,&i,1,&j,&scalar,INSERT_VALUES);CHKERRQ(ierr);</p>
<p>}</p>
<p>}</p>
<p>ierr=MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); </p>
<p>ierr=MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);</p>
<p>/* Inserting the non-zero elements. */</p>
<p>for(i=0;i<lenB;++i)</p>
<p>{</p>
<p>tempi=(int)B[i][0]-1;</p>
<p>tempj=(int)B[i][1]-1;</p>
<p>ierr=MatSetValues(A,1,&tempi,1,&tempj,&B[i][2],INSERT_VALUES);CHKERRQ(ierr);</p>
<p>}</p>
<p>ierr=MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); </p>
<p>ierr=MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);</p>
<p></p>
<p>//ierr=MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); </p>
<p>/*</p>
<p>Destroy any objects created.</p>
<p>*/</p>
<p>ierr=VecDestroy(b);CHKERRQ(ierr);</p>
<p></p>
<p>ierr=MatDestroy(A);CHKERRQ(ierr);</p>
<p>ierr=PetscFinalize();CHKERRQ(ierr);</p>
<p>return 0;</p>
<p>}</p></span></div>
<div> </div>
<div> </div>
<div> </div>
<div> </div>
<div> </div>
<div> </div>
<div> </div>
<div> </div>
<div> </div>
<div> </div>
<div> </div>
<div> </div>
</blockquote></div><br><br clear="all"><br>-- <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>
</div>