Hi,
I'm using (from error log) "Petsc Development HG revision:
0216b5820377f721a70342e04a44fa54f3883ec9 HG Date: Sun Oct 30 23:58:53 2011
-0500 "
I had some problems using MatMatMultTranspose for SeqAIJ matrices. It seems
to me that MatMatMultTranspose(A,B,...,MAT_INITIAL_MATRIX,&C), which should
perform C=A^T*B according to the web pages, doesn't its job. The
pseudo-code listed below prints
sizes A_II 64 64
sizes phi D 64 30
check sizes TM1 30 30
check sizes TM2 64 64
check sizes TM3 64 64
check sizes AUXMAT2 64 30
>From the code below, either TM2 and TM3 sizes are wrong. They should be
both 30 30. It seems that it flips matrices dimensions in creating the
output of MatMatMultTranspose.
An inspection of matrix values give me no further insights.
--------------------------------------------------------------------------------------------------------------------
Mat TM1,TM2,TM3,AUXMAT2;
Mat coarse_phi_D;
.....all matrices involved are Seq_AIJ matrices.............
ierr = MatGetSize(pcis->A_II,&rows,&columns);
printf("sizes A_II %d %d\n",rows,columns);
ierr = MatGetSize(coarse_phi_D,&rows,&columns);
printf("sizes phi D %d %d\n",rows,columns);
ierr =
MatPtAP(pcis->A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&TM1);
//TM1=PHI^T*A_II*PHI OK!
ierr =
MatMatMult(pcis->A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT2);
//AUXMAT2=A_II*PHI OK!
ierr =
MatMatMultTranspose(coarse_phi_D,AUXMAT2,MAT_INITIAL_MATRIX,1.0,&TM2);
//TM2=PHI^T*AUXMAT2 WRONG!
ierr =
MatMatMultTranspose(AUXMAT2,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&TM3);
//TM3=AUXMAT2^T*PHI WRONG!
ierr = MatGetSize(TM1,&rows,&columns);
printf("check sizes TM1 %d %d\n",rows,columns);
ierr = MatGetSize(TM2,&rows,&columns);
printf("check sizes TM2 %d %d\n",rows,columns);
ierr = MatGetSize(TM3,&rows,&columns);
printf("check sizes TM3 %d %d\n",rows,columns);
ierr = MatGetSize(AUXMAT2,&rows,&columns);
printf("check sizes AUXMAT2 %d %d\n",rows,columns);
----------------------------------------------------------------------------------------------------------------------------
--
Stefano
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
<http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20111109/aa7e1ed4/attachment.html>
