1 2 static char help[] = "Tests Elemental interface.\n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat Cdense,B,C,Ct; 9 Vec d; 10 PetscInt i,j,m = 5,n,nrows,ncols; 11 const PetscInt *rows,*cols; 12 IS isrows,iscols; 13 PetscErrorCode ierr; 14 PetscScalar *v; 15 PetscMPIInt rank,size; 16 PetscReal Cnorm; 17 PetscBool flg,mats_view=PETSC_FALSE; 18 19 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 20 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 21 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 22 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 23 n = m; 24 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 25 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view)); 26 27 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&C)); 28 CHKERRQ(MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE)); 29 CHKERRQ(MatSetType(C,MATELEMENTAL)); 30 CHKERRQ(MatSetFromOptions(C)); 31 CHKERRQ(MatSetUp(C)); 32 CHKERRQ(MatGetOwnershipIS(C,&isrows,&iscols)); 33 CHKERRQ(ISGetLocalSize(isrows,&nrows)); 34 CHKERRQ(ISGetIndices(isrows,&rows)); 35 CHKERRQ(ISGetLocalSize(iscols,&ncols)); 36 CHKERRQ(ISGetIndices(iscols,&cols)); 37 CHKERRQ(PetscMalloc1(nrows*ncols,&v)); 38 #if defined(PETSC_USE_COMPLEX) 39 PetscRandom rand; 40 PetscScalar rval; 41 CHKERRQ(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 42 CHKERRQ(PetscRandomSetFromOptions(rand)); 43 for (i=0; i<nrows; i++) { 44 for (j=0; j<ncols; j++) { 45 CHKERRQ(PetscRandomGetValue(rand,&rval)); 46 v[i*ncols+j] = rval; 47 } 48 } 49 CHKERRQ(PetscRandomDestroy(&rand)); 50 #else 51 for (i=0; i<nrows; i++) { 52 for (j=0; j<ncols; j++) { 53 v[i*ncols+j] = (PetscReal)(10000*rank+100*rows[i]+cols[j]); 54 } 55 } 56 #endif 57 CHKERRQ(MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES)); 58 CHKERRQ(ISRestoreIndices(isrows,&rows)); 59 CHKERRQ(ISRestoreIndices(iscols,&cols)); 60 CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 61 CHKERRQ(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 62 CHKERRQ(ISDestroy(&isrows)); 63 CHKERRQ(ISDestroy(&iscols)); 64 65 /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */ 66 CHKERRQ(MatDuplicate(C,MAT_COPY_VALUES,&B)); 67 if (mats_view) { 68 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Duplicated C:\n")); 69 CHKERRQ(MatView(B,PETSC_VIEWER_STDOUT_WORLD)); 70 } 71 CHKERRQ(MatDestroy(&B)); 72 CHKERRQ(MatConvert(C,MATMPIDENSE,MAT_INITIAL_MATRIX,&Cdense)); 73 CHKERRQ(MatMultEqual(C,Cdense,5,&flg)); 74 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Cdense != C. MatConvert() fails"); 75 76 /* Test MatNorm() */ 77 CHKERRQ(MatNorm(C,NORM_1,&Cnorm)); 78 79 /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */ 80 CHKERRQ(MatTranspose(C,MAT_INITIAL_MATRIX,&Ct)); 81 CHKERRQ(MatConjugate(Ct)); 82 if (mats_view) { 83 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"C's Transpose Conjugate:\n")); 84 CHKERRQ(MatView(Ct,PETSC_VIEWER_STDOUT_WORLD)); 85 } 86 CHKERRQ(MatZeroEntries(Ct)); 87 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&d)); 88 CHKERRQ(VecSetSizes(d,m>n ? n : m,PETSC_DECIDE)); 89 CHKERRQ(VecSetFromOptions(d)); 90 CHKERRQ(MatGetDiagonal(C,d)); 91 if (mats_view) { 92 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Diagonal of C:\n")); 93 CHKERRQ(VecView(d,PETSC_VIEWER_STDOUT_WORLD)); 94 } 95 if (m>n) { 96 CHKERRQ(MatDiagonalScale(C,NULL,d)); 97 } else { 98 CHKERRQ(MatDiagonalScale(C,d,NULL)); 99 } 100 if (mats_view) { 101 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Diagonal Scaled C:\n")); 102 CHKERRQ(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 103 } 104 105 /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */ 106 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&B)); 107 CHKERRQ(MatSetSizes(B,m,n,PETSC_DECIDE,PETSC_DECIDE)); 108 CHKERRQ(MatSetType(B,MATELEMENTAL)); 109 CHKERRQ(MatSetFromOptions(B)); 110 CHKERRQ(MatSetUp(B)); 111 CHKERRQ(MatGetOwnershipIS(B,&isrows,&iscols)); 112 CHKERRQ(ISGetLocalSize(isrows,&nrows)); 113 CHKERRQ(ISGetIndices(isrows,&rows)); 114 CHKERRQ(ISGetLocalSize(iscols,&ncols)); 115 CHKERRQ(ISGetIndices(iscols,&cols)); 116 for (i=0; i<nrows; i++) { 117 for (j=0; j<ncols; j++) { 118 v[i*ncols+j] = (PetscReal)(1000*rows[i]+cols[j]); 119 } 120 } 121 CHKERRQ(MatSetValues(B,nrows,rows,ncols,cols,v,INSERT_VALUES)); 122 CHKERRQ(PetscFree(v)); 123 CHKERRQ(ISRestoreIndices(isrows,&rows)); 124 CHKERRQ(ISRestoreIndices(iscols,&cols)); 125 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 126 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 127 CHKERRQ(MatAXPY(B,2.5,C,SAME_NONZERO_PATTERN)); 128 CHKERRQ(MatAYPX(B,3.75,C,SAME_NONZERO_PATTERN)); 129 CHKERRQ(MatConvert(B,MATDENSE,MAT_INPLACE_MATRIX,&B)); 130 if (mats_view) { 131 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"B after MatAXPY and MatAYPX:\n")); 132 CHKERRQ(MatView(B,PETSC_VIEWER_STDOUT_WORLD)); 133 } 134 CHKERRQ(ISDestroy(&isrows)); 135 CHKERRQ(ISDestroy(&iscols)); 136 CHKERRQ(MatDestroy(&B)); 137 138 /* Test MatMatTransposeMult(): B = C*C^T */ 139 CHKERRQ(MatMatTransposeMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B)); 140 CHKERRQ(MatScale(C,2.0)); 141 CHKERRQ(MatMatTransposeMult(C,C,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B)); 142 CHKERRQ(MatMatTransposeMultEqual(C,C,B,10,&flg)); 143 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatTransposeMult: B != C*B^T"); 144 145 if (mats_view) { 146 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"C MatMatTransposeMult C:\n")); 147 CHKERRQ(MatView(B,PETSC_VIEWER_STDOUT_WORLD)); 148 } 149 150 CHKERRQ(MatDestroy(&Cdense)); 151 CHKERRQ(PetscFree(v)); 152 CHKERRQ(MatDestroy(&B)); 153 CHKERRQ(MatDestroy(&C)); 154 CHKERRQ(MatDestroy(&Ct)); 155 CHKERRQ(VecDestroy(&d)); 156 ierr = PetscFinalize(); 157 return ierr; 158 } 159 160 /*TEST 161 162 test: 163 nsize: 2 164 args: -m 3 -n 2 165 requires: elemental 166 167 test: 168 suffix: 2 169 nsize: 6 170 args: -m 2 -n 3 171 requires: elemental 172 173 test: 174 suffix: 3 175 nsize: 1 176 args: -m 2 -n 3 177 requires: elemental 178 output_file: output/ex39_1.out 179 180 TEST*/ 181