I managed to do it by looking at: *PETScWrappers::MatrixBase*, but my
solution is extremely ugly!
double max_coup_dofs = dof_handler.max_couplings_between_dofs();
PETScWrappers::SparseMatrix local_M;
local_M.reinit (par.dofs,par.dofs, max_coup_dofs);
MatSetOption( local_M, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
{ // reading and parsing from (Mat M_read) into
PETScWrappers::SparseMatrix local_M;
PetscInt begin, end;
const int ierr = MatGetOwnershipRange (
static_cast<const Mat &>(M_read), &begin, &end);
std::pair<PETScWrappers::MatrixBase::size_type,
PETScWrappers::MatrixBase::size_type>
loc_range = std::make_pair (begin, end);
PetscInt ncols;
const PetscInt *colnums;
const PetscScalar *values;
PETScWrappers::MatrixBase::size_type row;
for (row = loc_range.first; row < loc_range.
second; ++row) {
int ierr = MatGetRow(M_read, row, &ncols, &
colnums, &values);
(void)ierr;
AssertThrow (ierr == 0, PETScWrappers::
MatrixBase::ExcPETScError(ierr));
for (PetscInt col = 0; col < ncols; ++col)
{
// cout << "(" << row << "," <<
colnums[col] << ") " << values[col] << std::endl;
local_M.set (row,colnums[col],
values[col] );
}
ierr = MatRestoreRow(M_read, row, &ncols, &
colnums, &values);
AssertThrow (ierr == 0, PETScWrappers::
MatrixBase::ExcPETScError(ierr));
}
//AssertThrow (out, ExcIO());
}
local_M.compress (VectorOperation::add);
//printf("\n\n\n ******** local_M norm is %f ******* \n\n\n",
local_M.l1_norm () );
if (abs (result-local_M.l1_norm ()) > 1e-10) {
printf("\n\n ******** M norm difference greater than 1e-10 after
reading ******* \n\n" );
} else {
printf("\n\n Mass matrix successfully read! \n\n" );
}
Any simpler alternative?
By the way, I noticed also the definition:
MatrixBase::operator Mat () const
{
return matrix;
}
By I was unable to use it!
Thanks,
Juan Carlos Araújo,
Umeå Universitet
El miércoles, 24 de mayo de 2017, 11:18:46 (UTC+2), Juan Carlos Araujo
Cabarcas escribió:
>
> Dear all,
>
> I have a matrix written by using PETSc in a Binary file, and I am trying
> work with it in my dealii environment. For this, I read a system matrix
> from PETSc:
>
> Mat M_read;
>
> PetscViewerBinaryOpen(PETSC_COMM_WORLD,file.c_str (),FILE_MODE_READ,&fd);
> MatCreate(PETSC_COMM_WORLD, &M_read);
> MatSetOptionsPrefix( M_read,"m_");
> MatSetFromOptions( M_read );
> MatLoad( local_M ,fd);
> PetscViewerDestroy(&fd);
> MatGetSize( M_read ,&mm,&nn);
> MatGetInfo( M_read ,MAT_LOCAL,&matinfo);
> printf("matinfo.nz_used %g\n", matinfo.nz_used);
>
> and I would like to initialize a *dealii PETScWrappers::SparseMatrix *by
> using *M_read*.
>
> If I try the naĩve way:
> PETScWrappers::SparseMatrix M = PETScWrappers::SparseMatrix (M_read);
>
> I get the error:
>
> error: ‘dealii::PETScWrappers::SparseMatrix::SparseMatrix(const
> dealii::PETScWrappers::SparseMatrix&)’ is private
> SparseMatrix(const SparseMatrix &);
>
> I know the *PETScWrappers::SparseMatrix* contains the protected member *Mat
> matrix*, but I quite don't know how to copy it to initialize *M*.
>
> Is there a way to achieve this? any help is appreciated!
>
> Thanks in advance,
> Juan Carlos Araújo,
> Umeå Universitet
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
