[petsc-users] Problem with MATSOLVERPASTIX using KSPSolveTranspose

[Zin Lin]

Hi,

I am having trouble calling KSPSolveTranspose under MATSOLVERPASTIX.
The following is a very simple C code solving A^T x = b to illustrate the
issue. The matrix A and vector b are built following the PETSC tutorial.
KSPSolve works fine but when it comes to KSPSolveTranspose, it spits out
the error:

[0]PETSC ERROR: --------------------- Error Message
------------------------------------
[0]PETSC ERROR: No support for this operation for this object type!
[0]PETSC ERROR: Matrix type mpiaij!

The code runs fine if I change the solver package MATSOLVERPASTIX to
MATSOLVERMUMPS.
However, I need the PASTIX solver and the transpose solution for solving a
large problem.
Could you please let me know what is wrong with the PASTIX solver or with
my usage?

Best regards,
Zin


/******************/
#include <stdlib.h>
#include <petsc.h>
#include <string.h>
#include <complex.h>

int main(int argc, char **argv)
{
 /* -------Initialize ------*/
PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
PetscPrintf(PETSC_COMM_WORLD,"--------Initializing------ \n");
PetscErrorCode ierr;

Vec x,y;
Mat A;
int N=10;
int low, high, i;
double val;
int start, end;
double v[3];
int row, cols[3];

ierr = VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE,N, &x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&y); CHKERRQ(ierr);

VecGetOwnershipRange(x, &low, &high);
val = low*10.0;
for(i = low; i < high; ++i) {
 VecSetValues(x, 1, &i, &val, INSERT_VALUES);
 val += 10.0;
}
VecAssemblyBegin(x);
VecAssemblyEnd(x);

MatCreate(PETSC_COMM_WORLD, &A);
MatSetType(A,MATMPIAIJ);
MatSetSizes(A,PETSC_DECIDE, PETSC_DECIDE, N, N);
MatMPIAIJSetPreallocation(A, 3, PETSC_NULL, 3, PETSC_NULL);
v[0] = -1.0; v[1] = 2.0; v[2] = -1.0;
MatGetOwnershipRange(A,&start,&end);
for(row = start; row < end; row++) {
 cols[0] = row-1; cols[1] = row; cols[2] = row+1;
if (row == 0) {
 MatSetValues(A,1,&row,2,&cols[1],&v[1],INSERT_VALUES);
} else if (row == N-1) {
 MatSetValues(A,1,&row,2,cols,v,INSERT_VALUES);
} else {
 MatSetValues(A,1,&row,3,cols,v,INSERT_VALUES);
}
}
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

PetscPrintf(PETSC_COMM_WORLD,"A and x setup. \n");

KSP ksp;
PC pc;
 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
 ierr = KSPSetType(ksp, KSPGMRES);CHKERRQ(ierr);
 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
 ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
 ierr = PCFactorSetMatSolverPackage(pc,MATSOLVERPASTIX);CHKERRQ(ierr);
 ierr = PCSetFromOptions(pc);
 ierr =
KSPSetTolerances(ksp,1e-14,PETSC_DEFAULT,PETSC_DEFAULT,20);CHKERRQ(ierr);
 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

 ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
 ierr = KSPSolve(ksp,x,y);CHKERRQ(ierr);

double p;
 ierr = VecDot(x,y,&p); CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"x dot y: %g. \n", p);

Vec u;
ierr = VecDuplicate(x,&u); CHKERRQ(ierr);
ierr = KSPSolveTranspose(ksp,y,u);CHKERRQ(ierr);

PetscPrintf(PETSC_COMM_WORLD,"Transpose solved. \n");

 ierr = VecDestroy(&x); CHKERRQ(ierr);
 ierr = VecDestroy(&y); CHKERRQ(ierr);
 ierr = MatDestroy(&A); CHKERRQ(ierr);
 ierr = KSPDestroy(&ksp); CHKERRQ(ierr);

 {
   int rank;
   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
   MPI_Barrier(PETSC_COMM_WORLD);
 }

 ierr = PetscFinalize(); CHKERRQ(ierr);

return 0;
}
/******************/
-- 
Zin Lin
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140725/981f1473/attachment.html>