1 static char help[] = "Test conversion of ScaLAPACK matrices.\n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char** argv) 6 { 7 Mat A,A_scalapack; 8 PetscInt i,j,M=10,N=5,nloc,mloc,nrows,ncols; 9 PetscErrorCode ierr; 10 PetscMPIInt rank,size; 11 IS isrows,iscols; 12 const PetscInt *rows,*cols; 13 PetscScalar *v; 14 MatType type; 15 PetscBool isDense,isAIJ,flg; 16 17 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 18 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 19 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 20 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 21 ierr = PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);CHKERRQ(ierr); 22 23 /* Create a matrix */ 24 ierr = MatCreate(PETSC_COMM_WORLD, &A);CHKERRQ(ierr); 25 mloc = PETSC_DECIDE; 26 ierr = PetscSplitOwnershipEqual(PETSC_COMM_WORLD,&mloc,&M);CHKERRQ(ierr); 27 nloc = PETSC_DECIDE; 28 ierr = PetscSplitOwnershipEqual(PETSC_COMM_WORLD,&nloc,&N);CHKERRQ(ierr); 29 ierr = MatSetSizes(A,mloc,nloc,M,N);CHKERRQ(ierr); 30 ierr = MatSetType(A,MATDENSE);CHKERRQ(ierr); 31 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 32 ierr = MatSetUp(A);CHKERRQ(ierr); 33 34 /* Set local matrix entries */ 35 ierr = MatGetOwnershipIS(A,&isrows,&iscols);CHKERRQ(ierr); 36 ierr = ISGetLocalSize(isrows,&nrows);CHKERRQ(ierr); 37 ierr = ISGetIndices(isrows,&rows);CHKERRQ(ierr); 38 ierr = ISGetLocalSize(iscols,&ncols);CHKERRQ(ierr); 39 ierr = ISGetIndices(iscols,&cols);CHKERRQ(ierr); 40 ierr = PetscMalloc1(nrows*ncols,&v);CHKERRQ(ierr); 41 42 for (i=0; i<nrows; i++) { 43 for (j=0; j<ncols; j++) { 44 if (size == 1) { 45 v[i*ncols+j] = (PetscScalar)(i+j); 46 } else { 47 v[i*ncols+j] = (PetscScalar)rank+j*0.1; 48 } 49 } 50 } 51 ierr = MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES);CHKERRQ(ierr); 52 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 53 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 54 55 /* Test MatSetValues() by converting A to A_scalapack */ 56 ierr = MatGetType(A,&type);CHKERRQ(ierr); 57 if (size == 1) { 58 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&isDense);CHKERRQ(ierr); 59 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isAIJ);CHKERRQ(ierr); 60 } else { 61 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIDENSE,&isDense);CHKERRQ(ierr); 62 ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isAIJ);CHKERRQ(ierr); 63 } 64 65 if (isDense || isAIJ) { 66 Mat Aexplicit; 67 ierr = MatConvert(A,MATSCALAPACK,MAT_INITIAL_MATRIX,&A_scalapack);CHKERRQ(ierr); 68 ierr = MatComputeOperator(A_scalapack,isAIJ?MATAIJ:MATDENSE,&Aexplicit);CHKERRQ(ierr); 69 ierr = MatMultEqual(Aexplicit,A_scalapack,5,&flg);CHKERRQ(ierr); 70 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Aexplicit != A_scalapack."); 71 ierr = MatDestroy(&Aexplicit);CHKERRQ(ierr); 72 73 /* Test MAT_REUSE_MATRIX which is only supported for inplace conversion */ 74 ierr = MatConvert(A,MATSCALAPACK,MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr); 75 ierr = MatMultEqual(A_scalapack,A,5,&flg);CHKERRQ(ierr); 76 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"A_scalapack != A."); 77 ierr = MatDestroy(&A_scalapack);CHKERRQ(ierr); 78 } 79 80 ierr = ISRestoreIndices(isrows,&rows);CHKERRQ(ierr); 81 ierr = ISRestoreIndices(iscols,&cols);CHKERRQ(ierr); 82 ierr = ISDestroy(&isrows);CHKERRQ(ierr); 83 ierr = ISDestroy(&iscols);CHKERRQ(ierr); 84 ierr = PetscFree(v);CHKERRQ(ierr); 85 ierr = MatDestroy(&A);CHKERRQ(ierr); 86 ierr = PetscFinalize(); 87 return ierr; 88 } 89 90 /*TEST 91 92 build: 93 requires: scalapack 94 95 test: 96 nsize: 6 97 98 test: 99 suffix: 2 100 nsize: 6 101 args: -mat_type aij 102 output_file: output/ex243_1.out 103 104 test: 105 suffix: 3 106 nsize: 6 107 args: -mat_type scalapack 108 output_file: output/ex243_1.out 109 110 TEST*/ 111