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