Dear all, I was trying to write a simple class template code to interface Petsc matrices along with Boost sparse matrices and vice versa. However, my first naive try did not work out of the box, since the template class requires some functions from PETSc libraries and these functions should be located in the object files to be able to compile. Since class template is header only and there is no object code generation, it can not find the PETSc library functions which is logical. I am also new to templates in C++. I just gave a try to create a template class instead of code duplication however it did not end up as expected. Are there some smarter ways to accomplish this task. The function analyzes a boost csr matrix and sets the rows of a Petsc matrix. T is for boost matrices and T1 for Mat in Petsc. I can also get some comments on the code ;) Best regards, Umut template <class T, class T1> int converter<T, T1>::convertMBo2Pe( const T& A, T1 A_ ){ PetscErrorCode ierr; int cntNnz = 0; typedef typename T::iterator1 i1_t; typedef typename T::iterator2 i2_t; //int nnz[ooFelieMatrix.size1()]; int nnz[A.size1()]; unsigned ind=0; //get information about the matrix double* vals = NULL; for (i1_t i1 = A.begin1(); i1 != A.end1(); ++i1) { nnz[ind] = distance(i1.begin(), i1.end()); ind++; } // create the matrix depending // on the values of the nonzeros // on each row ierr = MatCreateSeqAIJ( PETSC_COMM_SELF, A.size1(), A.size2(), cntNnz, nnz, A_ ); PetscInt rInd = 0, cInd=0; PetscInt* rCount, dummy; rCount = &dummy; // pointer to values in a row PetscScalar* valsOfRowI = NULL; PetscInt* colIndexOfRowI = NULL; PetscInt rC = 1; for(i1_t i1 = A.begin1(); i1 != A.end1(); ++i1) { // allocate space for the values of row I valsOfRowI = new PetscScalar[nnz[rInd]]; colIndexOfRowI = new PetscInt[nnz[rInd]]; for(i2_t i2 = i1.begin(); i2 != i1.end(); ++i2) { colIndexOfRowI[cInd] = i2.index2(); valsOfRowI[cInd] = *i2; cInd++; } // setting one row each time *rCount = rInd; MatSetValues( A_, rC, rCount, nnz[rInd], colIndexOfRowI, valsOfRowI, INSERT_VALUES ); // delete delete [] valsOfRowI; delete [] colIndexOfRowI; rInd++; cInd = 0; } // MatAssemblyBegin( A_, MAT_FINAL_ASSEMBLY ); MatAssemblyEnd( A_, MAT_FINAL_ASSEMBLY ); // return return 0; } -- - Hope is a good thing, maybe the best of things and no good thing ever dies... The Shawshank Redemption, replique of Tim Robbins
What is the error? Matt On Mon, Dec 13, 2010 at 3:11 PM, Umut Tabak <[email protected]> wrote:
Dear all,
I was trying to write a simple class template code to interface Petsc matrices along with Boost sparse matrices and vice versa. However, my first naive try did not work out of the box, since the template class requires some functions from PETSc libraries and these functions should be located in the object files to be able to compile. Since class template is header only and there is no object code generation, it can not find the PETSc library functions which is logical. I am also new to templates in C++. I just gave a try to create a template class instead of code duplication however it did not end up as expected. Are there some smarter ways to accomplish this task. The function analyzes a boost csr matrix and sets the rows of a Petsc matrix. T is for boost matrices and T1 for Mat in Petsc. I can also get some comments on the code ;)
Best regards, Umut
template <class T, class T1> int converter<T, T1>::convertMBo2Pe( const T& A, T1 A_ ){ PetscErrorCode ierr; int cntNnz = 0; typedef typename T::iterator1 i1_t; typedef typename T::iterator2 i2_t; //int nnz[ooFelieMatrix.size1()]; int nnz[A.size1()]; unsigned ind=0; //get information about the matrix
double* vals = NULL; for (i1_t i1 = A.begin1(); i1 != A.end1(); ++i1) { nnz[ind] = distance(i1.begin(), i1.end()); ind++; } // create the matrix depending // on the values of the nonzeros // on each row ierr = MatCreateSeqAIJ( PETSC_COMM_SELF, A.size1(), A.size2(), cntNnz, nnz, A_ ); PetscInt rInd = 0, cInd=0; PetscInt* rCount, dummy; rCount = &dummy; // pointer to values in a row PetscScalar* valsOfRowI = NULL; PetscInt* colIndexOfRowI = NULL; PetscInt rC = 1; for(i1_t i1 = A.begin1(); i1 != A.end1(); ++i1) { // allocate space for the values of row I valsOfRowI = new PetscScalar[nnz[rInd]]; colIndexOfRowI = new PetscInt[nnz[rInd]]; for(i2_t i2 = i1.begin(); i2 != i1.end(); ++i2) { colIndexOfRowI[cInd] = i2.index2(); valsOfRowI[cInd] = *i2; cInd++; } // setting one row each time *rCount = rInd; MatSetValues( A_, rC, rCount, nnz[rInd], colIndexOfRowI, valsOfRowI, INSERT_VALUES ); // delete delete [] valsOfRowI; delete [] colIndexOfRowI; rInd++; cInd = 0; } // MatAssemblyBegin( A_, MAT_FINAL_ASSEMBLY ); MatAssemblyEnd( A_, MAT_FINAL_ASSEMBLY ); // return return 0; }
-- - Hope is a good thing, maybe the best of things and no good thing ever dies... The Shawshank Redemption, replique of Tim Robbins
-- 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
participants (2)
-
Matthew Knepley -
Umut Tabak