1 2 static char help[] = "Test interface of Elemental. \n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat C,Caij; 9 PetscInt i,j,m = 5,n,nrows,ncols; 10 const PetscInt *rows,*cols; 11 IS isrows,iscols; 12 PetscErrorCode ierr; 13 PetscBool flg,Test_MatMatMult=PETSC_FALSE,mats_view=PETSC_FALSE; 14 PetscScalar *v; 15 PetscMPIInt rank,size; 16 17 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 18 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 19 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 20 ierr = PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view);CHKERRQ(ierr); 21 22 /* Get local block or element size*/ 23 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 24 n = m; 25 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 26 27 ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr); 28 ierr = MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 29 ierr = MatSetType(C,MATELEMENTAL);CHKERRQ(ierr); 30 ierr = MatSetFromOptions(C);CHKERRQ(ierr); 31 ierr = MatSetUp(C);CHKERRQ(ierr); 32 33 ierr = PetscOptionsHasName(NULL,NULL,"-row_oriented",&flg);CHKERRQ(ierr); 34 if (flg) {ierr = MatSetOption(C,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr);} 35 ierr = MatGetOwnershipIS(C,&isrows,&iscols);CHKERRQ(ierr); 36 ierr = PetscOptionsHasName(NULL,NULL,"-Cexp_view_ownership",&flg);CHKERRQ(ierr); 37 if (flg) { /* View ownership of explicit C */ 38 IS tmp; 39 ierr = PetscPrintf(PETSC_COMM_WORLD,"Ownership of explicit C:\n");CHKERRQ(ierr); 40 ierr = PetscPrintf(PETSC_COMM_WORLD,"Row index set:\n");CHKERRQ(ierr); 41 ierr = ISOnComm(isrows,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp);CHKERRQ(ierr); 42 ierr = ISView(tmp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 43 ierr = ISDestroy(&tmp);CHKERRQ(ierr); 44 ierr = PetscPrintf(PETSC_COMM_WORLD,"Column index set:\n");CHKERRQ(ierr); 45 ierr = ISOnComm(iscols,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp);CHKERRQ(ierr); 46 ierr = ISView(tmp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 47 ierr = ISDestroy(&tmp);CHKERRQ(ierr); 48 } 49 50 /* Set local matrix entries */ 51 ierr = ISGetLocalSize(isrows,&nrows);CHKERRQ(ierr); 52 ierr = ISGetIndices(isrows,&rows);CHKERRQ(ierr); 53 ierr = ISGetLocalSize(iscols,&ncols);CHKERRQ(ierr); 54 ierr = ISGetIndices(iscols,&cols);CHKERRQ(ierr); 55 ierr = PetscMalloc1(nrows*ncols,&v);CHKERRQ(ierr); 56 for (i=0; i<nrows; i++) { 57 for (j=0; j<ncols; j++) { 58 /*v[i*ncols+j] = (PetscReal)(rank);*/ 59 v[i*ncols+j] = (PetscReal)(rank*10000+100*rows[i]+cols[j]); 60 } 61 } 62 ierr = MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES);CHKERRQ(ierr); 63 ierr = ISRestoreIndices(isrows,&rows);CHKERRQ(ierr); 64 ierr = ISRestoreIndices(iscols,&cols);CHKERRQ(ierr); 65 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 66 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 67 68 /* Test MatView() */ 69 if (mats_view) { 70 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 71 } 72 73 /* Test MatMissingDiagonal() */ 74 ierr = MatMissingDiagonal(C,&flg,NULL);CHKERRQ(ierr); 75 if (flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"MatMissingDiagonal() did not return false!\n"); 76 77 /* Set unowned matrix entries - add subdiagonals and diagonals from proc[0] */ 78 if (!rank) { 79 PetscInt M,N,cols[2]; 80 ierr = MatGetSize(C,&M,&N);CHKERRQ(ierr); 81 for (i=0; i<M; i++) { 82 cols[0] = i; v[0] = i + 0.5; 83 cols[1] = i-1; v[1] = 0.5; 84 if (i) { 85 ierr = MatSetValues(C,1,&i,2,cols,v,ADD_VALUES);CHKERRQ(ierr); 86 } else { 87 ierr = MatSetValues(C,1,&i,1,&i,v,ADD_VALUES);CHKERRQ(ierr); 88 } 89 } 90 } 91 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 92 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 93 94 /* Test MatMult() */ 95 ierr = MatComputeOperator(C,MATAIJ,&Caij);CHKERRQ(ierr); 96 ierr = MatMultEqual(C,Caij,5,&flg);CHKERRQ(ierr); 97 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultEqual() fails"); 98 ierr = MatMultTransposeEqual(C,Caij,5,&flg);CHKERRQ(ierr); 99 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeEqual() fails"); 100 101 /* Test MatMultAdd() and MatMultTransposeAddEqual() */ 102 ierr = MatMultAddEqual(C,Caij,5,&flg);CHKERRQ(ierr); 103 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultAddEqual() fails"); 104 ierr = MatMultTransposeAddEqual(C,Caij,5,&flg);CHKERRQ(ierr); 105 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeAddEqual() fails"); 106 107 /* Test MatMatMult() */ 108 ierr = PetscOptionsHasName(NULL,NULL,"-test_matmatmult",&Test_MatMatMult);CHKERRQ(ierr); 109 if (Test_MatMatMult) { 110 Mat CCelem,CCaij; 111 ierr = MatMatMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCelem);CHKERRQ(ierr); 112 ierr = MatMatMult(Caij,Caij,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCaij);CHKERRQ(ierr); 113 ierr = MatMultEqual(CCelem,CCaij,5,&flg);CHKERRQ(ierr); 114 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"CCelem != CCaij. MatMatMult() fails"); 115 ierr = MatDestroy(&CCaij);CHKERRQ(ierr); 116 ierr = MatDestroy(&CCelem);CHKERRQ(ierr); 117 } 118 119 ierr = PetscFree(v);CHKERRQ(ierr); 120 ierr = ISDestroy(&isrows);CHKERRQ(ierr); 121 ierr = ISDestroy(&iscols);CHKERRQ(ierr); 122 ierr = MatDestroy(&C);CHKERRQ(ierr); 123 ierr = MatDestroy(&Caij);CHKERRQ(ierr); 124 ierr = PetscFinalize(); 125 return ierr; 126 } 127 128 129 130 131 /*TEST 132 133 test: 134 nsize: 2 135 requires: elemental 136 args: -mat_type elemental -m 2 -n 3 137 138 test: 139 suffix: 2 140 nsize: 6 141 requires: elemental 142 args: -mat_type elemental -m 2 -n 2 143 144 test: 145 suffix: 3 146 nsize: 6 147 requires: elemental 148 args: -mat_type elemental -m 2 -n 2 -test_matmatmult 149 output_file: output/ex38_2.out 150 151 TEST*/ 152