[petsc-users] Matrix Construction Question

[Jed Brown]

On Wed, Jul 6, 2011 at 09:53, Adam Byrd <adam1.byrd at gmail.com> wrote:

> I do not want the transpose, but that is what is coming out as the
> solution. Both my tests show the correct original matrix, and both yield the
> transpose of the solution verified against Mathematica and the website from
> which I pulled the 4x4 matrix in test.cpp. Could it be related to the
> factoring?


The problem is that you are trying to solve with an AIJ matrix instead of a
DENSE matrix. MatMatSolve() calls MatGetArray(). Since this is serial and
you preallocated enough for a dense matrix, it goes through with no memory
errors. The storage format for MATAIJ (that happens to be dense) and
MATDENSE is transposed, so you see the transpose. If you make the solution
matrix dense, then it will work correctly.


--- a/test.cpp
+++ b/test.cpp
@@ -48,7 +48,9 @@ int main(int argc,char **args)
        ierr = MatAssemblyEnd(identityMat,
MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

        //Setup inverseMat
-       ierr = MatDuplicate(testMat, MAT_DO_NOT_COPY_VALUES,
&inverseMat);CHKERRQ(ierr);
+       ierr = MatCreate(PETSC_COMM_WORLD, &inverseMat);CHKERRQ(ierr);
+       ierr = MatSetSizes(inverseMat, PETSC_DECIDE, PETSC_DECIDE, rows,
columns);CHKERRQ(ierr);
+       ierr = MatSetType(inverseMat, MATDENSE);CHKERRQ(ierr);
        ierr = MatAssemblyBegin(inverseMat,
MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
        ierr = MatAssemblyEnd(inverseMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);



PETSc should be more careful about matching types here.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110706/c480e964/attachment.htm>