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 <[email protected]>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