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 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 21 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 22 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 23 n = m; 24 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 25 ierr = PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view);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 ierr = MatGetOwnershipIS(C,&isrows,&iscols);CHKERRQ(ierr); 33 ierr = ISGetLocalSize(isrows,&nrows);CHKERRQ(ierr); 34 ierr = ISGetIndices(isrows,&rows);CHKERRQ(ierr); 35 ierr = ISGetLocalSize(iscols,&ncols);CHKERRQ(ierr); 36 ierr = ISGetIndices(iscols,&cols);CHKERRQ(ierr); 37 ierr = PetscMalloc1(nrows*ncols,&v);CHKERRQ(ierr); 38 #if defined(PETSC_USE_COMPLEX) 39 PetscRandom rand; 40 PetscScalar rval; 41 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr); 42 ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr); 43 for (i=0; i<nrows; i++) { 44 for (j=0; j<ncols; j++) { 45 ierr = PetscRandomGetValue(rand,&rval);CHKERRQ(ierr); 46 v[i*ncols+j] = rval; 47 } 48 } 49 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 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 ierr = MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES);CHKERRQ(ierr); 58 ierr = ISRestoreIndices(isrows,&rows);CHKERRQ(ierr); 59 ierr = ISRestoreIndices(iscols,&cols);CHKERRQ(ierr); 60 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 61 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 62 ierr = ISDestroy(&isrows);CHKERRQ(ierr); 63 ierr = ISDestroy(&iscols);CHKERRQ(ierr); 64 65 /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */ 66 ierr = MatDuplicate(C,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 67 if (mats_view) { 68 ierr = PetscPrintf(PETSC_COMM_WORLD,"Duplicated C:\n");CHKERRQ(ierr); 69 ierr = MatView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 70 } 71 ierr = MatDestroy(&B);CHKERRQ(ierr); 72 ierr = MatConvert(C,MATMPIDENSE,MAT_INITIAL_MATRIX,&Cdense);CHKERRQ(ierr); 73 ierr = MatMultEqual(C,Cdense,5,&flg);CHKERRQ(ierr); 74 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Cdense != C. MatConvert() fails"); 75 76 /* Test MatNorm() */ 77 ierr = MatNorm(C,NORM_1,&Cnorm);CHKERRQ(ierr); 78 79 /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */ 80 ierr = MatTranspose(C,MAT_INITIAL_MATRIX,&Ct);CHKERRQ(ierr); 81 ierr = MatConjugate(Ct);CHKERRQ(ierr); 82 if (mats_view) { 83 ierr = PetscPrintf(PETSC_COMM_WORLD,"C's Transpose Conjugate:\n");CHKERRQ(ierr); 84 ierr = MatView(Ct,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 85 } 86 ierr = MatZeroEntries(Ct);CHKERRQ(ierr); 87 ierr = VecCreate(PETSC_COMM_WORLD,&d);CHKERRQ(ierr); 88 ierr = VecSetSizes(d,m>n ? n : m,PETSC_DECIDE);CHKERRQ(ierr); 89 ierr = VecSetFromOptions(d);CHKERRQ(ierr); 90 ierr = MatGetDiagonal(C,d);CHKERRQ(ierr); 91 if (mats_view) { 92 ierr = PetscPrintf(PETSC_COMM_WORLD,"Diagonal of C:\n");CHKERRQ(ierr); 93 ierr = VecView(d,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 94 } 95 if (m>n) { 96 ierr = MatDiagonalScale(C,NULL,d);CHKERRQ(ierr); 97 } else { 98 ierr = MatDiagonalScale(C,d,NULL);CHKERRQ(ierr); 99 } 100 if (mats_view) { 101 ierr = PetscPrintf(PETSC_COMM_WORLD,"Diagonal Scaled C:\n");CHKERRQ(ierr); 102 ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 103 } 104 105 /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */ 106 ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr); 107 ierr = MatSetSizes(B,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 108 ierr = MatSetType(B,MATELEMENTAL);CHKERRQ(ierr); 109 ierr = MatSetFromOptions(B);CHKERRQ(ierr); 110 ierr = MatSetUp(B);CHKERRQ(ierr); 111 ierr = MatGetOwnershipIS(B,&isrows,&iscols);CHKERRQ(ierr); 112 ierr = ISGetLocalSize(isrows,&nrows);CHKERRQ(ierr); 113 ierr = ISGetIndices(isrows,&rows);CHKERRQ(ierr); 114 ierr = ISGetLocalSize(iscols,&ncols);CHKERRQ(ierr); 115 ierr = ISGetIndices(iscols,&cols);CHKERRQ(ierr); 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 ierr = MatSetValues(B,nrows,rows,ncols,cols,v,INSERT_VALUES);CHKERRQ(ierr); 122 ierr = PetscFree(v);CHKERRQ(ierr); 123 ierr = ISRestoreIndices(isrows,&rows);CHKERRQ(ierr); 124 ierr = ISRestoreIndices(iscols,&cols);CHKERRQ(ierr); 125 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 126 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 127 ierr = MatAXPY(B,2.5,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 128 ierr = MatAYPX(B,3.75,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 129 ierr = MatConvert(B,MATDENSE,MAT_INPLACE_MATRIX,&B);CHKERRQ(ierr); 130 if (mats_view) { 131 ierr = PetscPrintf(PETSC_COMM_WORLD,"B after MatAXPY and MatAYPX:\n");CHKERRQ(ierr); 132 ierr = MatView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 133 } 134 ierr = ISDestroy(&isrows);CHKERRQ(ierr); 135 ierr = ISDestroy(&iscols);CHKERRQ(ierr); 136 ierr = MatDestroy(&B);CHKERRQ(ierr); 137 138 /* Test MatMatTransposeMult(): B = C*C^T */ 139 ierr = MatMatTransposeMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 140 ierr = MatScale(C,2.0);CHKERRQ(ierr); 141 ierr = MatMatTransposeMult(C,C,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 142 ierr = MatMatTransposeMultEqual(C,C,B,10,&flg);CHKERRQ(ierr); 143 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatTransposeMult: B != C*B^T"); 144 145 if (mats_view) { 146 ierr = PetscPrintf(PETSC_COMM_WORLD,"C MatMatTransposeMult C:\n");CHKERRQ(ierr); 147 ierr = MatView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 148 } 149 150 ierr = MatDestroy(&Cdense);CHKERRQ(ierr); 151 ierr = PetscFree(v);CHKERRQ(ierr); 152 ierr = MatDestroy(&B);CHKERRQ(ierr); 153 ierr = MatDestroy(&C);CHKERRQ(ierr); 154 ierr = MatDestroy(&Ct);CHKERRQ(ierr); 155 ierr = VecDestroy(&d);CHKERRQ(ierr); 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