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