Re: [Boost-users] Parsing Sparse Matrix to Get JA, IA, A Vector

Subject: Re: [Boost-users] Parsing Sparse Matrix to Get JA, IA, A Vector
From: Damien Hocking (damien_at_[hidden])
Date: 2009-08-06 22:10:37

Peverall Dubois wrote:
> Does boost has a library to parse a sparse matrix
> to get the index the non-zero element value?
>
> E.g. given this matrix:
>
> 1 2 0 0
> 0 3 9 0
> 0 1 4 0
>
> We hope it to return
> A  = [ 1 2 3 9 1 4 ]
> IA = [ 1 3 5 7 ]
> JA = [ 1 2 2 3 2 3 ]
>
> Let NNZ denote the number of nonzero entries of Matrix. The first
> array is A, which is of length NNZ, and holds all nonzero entries of M
> in left-to-right top-to-bottom order. The second array is IA, which is
> of length m + 1 (i.e., one entry per row, plus one). IA(i) contains
> the index in A of the first nonzero element of row i. Row i of the
> original matrix extends from A(IA(i)) to A(IA(i+1)-1). The third
> array, JA, contains the column index of each element of A, so it also
> is of length NNZ.
>
> Is there any implementation for handling such a large matrix to give
> us A, IA and JA, without having have to slurp all the files into RAM?
> Because the actual dimension is very large 10^6 to 10^7.
>
> Regards,
> G.V.
> _______________________________________________
> Boost-users mailing list
> Boost-users_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/boost-users
>   
Trying to avoid top-posting... :-)

That's standard sparse Harwell-Boeing (or RB) format.  Are you trying to 
read from a file?  I have a file reader that I use for matrices from the 
Florida Sparse matrix collection that reads those files line by line.  
It reads all the extra descriptive parts that give the dimensional 
information about the matrix.  The memory usage isn't too bad.  Here's 
the code if you want it (not Boost quality, but it works).  Just call 
Readmatrix from main.

Damien

#include <iostream>
#include <fstream>
#include <string>
#include <stdlib.h>
using namespace std;

void finish(ifstream& file, char* line1, char* line2, char* line3, char* 
line4, char* line5)
{
   file.close();
   delete[]line1;
   delete[]line2;
   delete[]line3;
   delete[]line4;
   delete[]line5;
}

int readline(ifstream& file, char* line, int numchars)
{
   if (line == 0)
   {
      return 0;
   }
   int ctr = 0;
   cout << "\nReading line:\n";
   while(ctr < numchars)
   {
      file.read(&line[ctr], 1);
      cout << line[ctr];
      if (line[ctr] == '\n' || line[ctr] == '\r')
      {
         cout << "Read a newline at " << ctr;
         break;
      }
      ctr++;
   }
   return ctr;
}

void printline(char* line, int numchars)
{
   cout << "\n";
   for (int i = 0; i < numchars+1; i++)
   {
      cout << line[i];
   }
}

bool ReadMatrix(char* matrixfile, int* nrow, int* ncol, int* nonz, 
double** nzval, int** rowind, int** colptr, int* nrhs, double** rhs)
{
   ifstream inputfile(matrixfile);
   char* line1 = new char[81];
   char* line2 = new char[71];
   char* line3 = new char[71];
   char* line4 = new char[73];
   char* line5 = new char[43];
   inputfile.seekg(0);//reset file pointer to start of file
   cout << "\n";
   int charsread = 0;
   charsread = readline(inputfile, line1, 81);
   cout << "\n Matrix title is:\n";
   printline(line1, charsread);
   charsread = readline(inputfile, line2, 71);
   char* totLinesStr = &line2[0];
   int totLines = atoi(totLinesStr);
   char* ptrLinesStr = &line2[14];
   int ptrLines = atoi(ptrLinesStr);
   char* varIndexLinesStr = &line2[28];
   int varIndexLines = atoi(varIndexLinesStr);
   char* valueLinesStr = &line2[43];
   int valueLines = atoi(valueLinesStr);
   char* rhsLinesStr = &line2[56];
   int rhsLines = atoi(rhsLinesStr);
   cout << "\n";
   cout << "\nThere are " << totLines << " total lines of data.";
   cout << "\nThere are " << ptrLines << " lines of column pointers.";
   cout << "\nThere are " << varIndexLines << " lines of variable index 
values.";
   cout << "\nThere are " << valueLines << " lines of matrix values.";
   if (rhsLines > 0)
   {
      cout << "\nThere are " << rhsLines << " lines of RHS values.";
   }
   else
   {
      cout << "\nThere are no RHS values, so this data is for 
factorisation only, with no solution.";
   }
   int sum = ptrLines + varIndexLines + valueLines + rhsLines;
   cout << "\n\nColumn pointer lines + variable index lines + matrix 
value lines + RHS lines ";
   /*
   if (sum == totLines)
   {
   cout << " = Total lines";
   cout << "\nData is consistent, I can continue...";
   }
   else
   {
   cout << " does not equal Total lines";
   cout << "\nData is INCONSISTENT.  I cannot process further.";
   cout << "\nYou are the weakest link.  Goodbye.";
   finish(inputfile, line1, line2, line3, line4, line5);
   }
   */
   charsread = readline(inputfile, line3, 71);
   printline(line3, charsread);
   char matrixType[4];
   matrixType[0] = line3[0];
   matrixType[1] = line3[1];
   matrixType[2] = line3[2];
   matrixType[3] = '\0';
   if (strcmp(matrixType, "rua") == 0 || strcmp(matrixType, "RUA") == 0)
   {
      cout << "\n Matrix type is " << matrixType << ", OK";
   }
   else
   {
      cout << "\n Matrix type is " << matrixType << ", NOT OK";
      finish(inputfile, line1, line2, line3, line4, line5);
   }
   char* numRowsStr = &line3[14];
   int numRows = atoi(numRowsStr);
   *nrow = numRows;
   char* numColsStr = &line3[28];
   int numCols = atoi(numColsStr);
   *ncol = numCols;
   char* numEntriesStr = &line3[42];
   int numEntries = atoi(numEntriesStr);
   *nonz = numEntries;
   char* numElemEntriesStr = &line3[56];
   int numElemEntries = atoi(numElemEntriesStr);
   cout << "\n\nThe matrix has " << numRows << " rows.";
   cout << "\nThe matrix has " << numCols << " columns.";
   cout << "\nThe matrix has " << numEntries << " nonzeroes.\n";
   int numRHSVecs = 0;
   int numRowIndicesDummy = 0;
   charsread = readline(inputfile, line4, 73);
   bool haveGuess = false;
   bool haveAnswer = false;
   if (rhsLines > 0)
   {
      charsread = readline(inputfile, line5, 43);
      char rhsFormat = line5[0];
      char startGuess = line5[1];
      char answer = line5[2];
      if (rhsFormat == 'F')
      {
         char* numRHSVecStr = &line5[14];
         numRHSVecs = atoi(numRHSVecStr);
         cout << "\nI have " << numRHSVecs << " RHS vector(s), this 
matrix can be solved.";
         if (startGuess == 'G')
         {
            cout << "\n Read a G, I have an initial guess for the 
solution.";
            haveGuess = true;
         }
         if (answer== 'X')
         {
            cout << "\nI have a solution vector, solution can be checked.";
            haveAnswer = true;
         }
      }
      else
      {
         cout << "\nCan't read M format";
         finish(inputfile, line1, line2, line3, line4, line5);
      }
   } 
  
   int* colvalues = *colptr;
   (*colptr) = new int[numCols + 1];
   cout << "\n Reading column data (colptr)...\n"; 
   int ctr = 0;
   for (int i = 0; i < numCols + 1; i++)
   {
      inputfile >> (*colptr)[i];
      std::cout << " " << (*colptr)[i];
      ctr = i;
   }
   cout << "\n Last col value was at " << ctr << " value was " << 
(*colptr)[ctr];
  
   int* rowvalues = *rowind;
   (*rowind) = new int[numEntries];
   cout << "\n Reading row data (rowind)...\n"; 
   ctr = 0;
   for (int i = 0; i < numEntries; i++)
   {
      inputfile >> (*rowind)[i];
      std::cout << " " << (*rowind)[i];
      ctr = i;
   }
   cout << "\n Last row value was at " << ctr << " value was " << 
(*rowind)[ctr];
  
   double* matvalues = *nzval;
   (*nzval) = new double[numEntries];
   cout << "\n Reading matrix values...\n"; 
   ctr = 0;
   for (int i = 0; i < numEntries; i++)
   {
      inputfile >> (*nzval)[i];
      std::cout << " " << (*nzval)[i];
      ctr = i;
   }
  
   int nrhsvals = 0;
   inputfile >> nrhsvals;
   cout << "\n Last matrix value was at " << ctr << " value was " << 
(*nzval)[ctr];
  
   return true;  
}

//this is main:

int main(int argc, char *argv[])
{
   double *a;
   int *ia, *ja;
   int nrhs= 1, m, n, nnz;
   double *rhs, *x;
   double *a_mm;
   int* row_mm;
   int* col_mm;
  
   ReadMatrix("lhr17c.rb", &m, &n, &nnz, &a, &ia, &ja, &nrhs, &rhs);
   return 1;
}