1 /*$Id: axpy.c,v 1.54 2001/08/06 21:16:10 bsmith Exp $*/ 2 3 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "MatAXPY" 7 /*@ 8 MatAXPY - Computes Y = a*X + Y. 9 10 Collective on Mat 11 12 Input Parameters: 13 + X, Y - the matrices 14 - a - the PetscScalar multiplier 15 16 Contributed by: Matthew Knepley 17 18 Notes: 19 Since the current implementation of MatAXPY() uses MatGetRow() to access 20 matrix data, efficiency is somewhat limited. 21 22 Level: intermediate 23 24 .keywords: matrix, add 25 26 .seealso: MatAYPX() 27 @*/ 28 int MatAXPY(PetscScalar *a,Mat X,Mat Y) 29 { 30 int m1,m2,n1,n2,i,*row,start,end,j,ncols,ierr; 31 PetscScalar *val,*vals; 32 33 PetscFunctionBegin; 34 PetscValidHeaderSpecific(X,MAT_COOKIE); 35 PetscValidHeaderSpecific(Y,MAT_COOKIE); 36 PetscValidScalarPointer(a); 37 38 ierr = MatGetSize(X,&m1,&n1);CHKERRQ(ierr); 39 ierr = MatGetSize(Y,&m2,&n2);CHKERRQ(ierr); 40 if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_ERR_ARG_SIZ,"Non conforming matrix add: %d %d %d %d",m1,m2,n1,n2); 41 42 if (X->ops->axpy) { 43 ierr = (*X->ops->axpy)(a,X,Y);CHKERRQ(ierr); 44 } else { 45 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 46 if (*a == 1.0) { 47 for (i = start; i < end; i++) { 48 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 49 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 50 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 51 } 52 } else { 53 ierr = PetscMalloc((n1+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 54 for (i=start; i<end; i++) { 55 ierr = MatGetRow(X,i,&ncols,&row,&val);CHKERRQ(ierr); 56 for (j=0; j<ncols; j++) { 57 vals[j] = (*a)*val[j]; 58 } 59 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 60 ierr = MatRestoreRow(X,i,&ncols,&row,&val);CHKERRQ(ierr); 61 } 62 ierr = PetscFree(vals);CHKERRQ(ierr); 63 } 64 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 65 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 66 } 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "MatShift" 72 /*@ 73 MatShift - Computes Y = Y + a I, where a is a PetscScalar and I is the identity matrix. 74 75 Collective on Mat 76 77 Input Parameters: 78 + Y - the matrices 79 - a - the PetscScalar 80 81 Level: intermediate 82 83 .keywords: matrix, add, shift 84 85 .seealso: MatDiagonalSet() 86 @*/ 87 int MatShift(PetscScalar *a,Mat Y) 88 { 89 int i,start,end,ierr; 90 91 PetscFunctionBegin; 92 PetscValidHeaderSpecific(Y,MAT_COOKIE); 93 PetscValidScalarPointer(a); 94 if (Y->ops->shift) { 95 ierr = (*Y->ops->shift)(a,Y);CHKERRQ(ierr); 96 } else { 97 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 98 for (i=start; i<end; i++) { 99 ierr = MatSetValues(Y,1,&i,1,&i,a,ADD_VALUES);CHKERRQ(ierr); 100 } 101 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 102 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 103 } 104 PetscFunctionReturn(0); 105 } 106 107 #undef __FUNCT__ 108 #define __FUNCT__ "MatDiagonalSet" 109 /*@ 110 MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix 111 that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is 112 INSERT_VALUES. 113 114 Input Parameters: 115 + Y - the input matrix 116 . D - the diagonal matrix, represented as a vector 117 - i - INSERT_VALUES or ADD_VALUES 118 119 Collective on Mat and Vec 120 121 Level: intermediate 122 123 .keywords: matrix, add, shift, diagonal 124 125 .seealso: MatShift() 126 @*/ 127 int MatDiagonalSet(Mat Y,Vec D,InsertMode is) 128 { 129 int i,start,end,ierr; 130 131 PetscFunctionBegin; 132 PetscValidHeaderSpecific(Y,MAT_COOKIE); 133 PetscValidHeaderSpecific(D,VEC_COOKIE); 134 if (Y->ops->diagonalset) { 135 ierr = (*Y->ops->diagonalset)(Y,D,is);CHKERRQ(ierr); 136 } else { 137 int vstart,vend; 138 PetscScalar *v; 139 ierr = VecGetOwnershipRange(D,&vstart,&vend);CHKERRQ(ierr); 140 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 141 if (vstart != start || vend != end) { 142 SETERRQ4(PETSC_ERR_ARG_SIZ,"Vector ownership range not compatible with matrix: %d %d vec %d %d mat",vstart,vend,start,end); 143 } 144 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 145 for (i=start; i<end; i++) { 146 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,is);CHKERRQ(ierr); 147 } 148 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 149 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 150 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 151 } 152 PetscFunctionReturn(0); 153 } 154 155 #undef __FUNCT__ 156 #define __FUNCT__ "MatAYPX" 157 /*@ 158 MatAYPX - Computes Y = X + a*Y. 159 160 Collective on Mat 161 162 Input Parameters: 163 + X,Y - the matrices 164 - a - the PetscScalar multiplier 165 166 Contributed by: Matthew Knepley 167 168 Notes: 169 This routine currently uses the MatAXPY() implementation. 170 171 Level: intermediate 172 173 .keywords: matrix, add 174 175 .seealso: MatAXPY() 176 @*/ 177 int MatAYPX(PetscScalar *a,Mat X,Mat Y) 178 { 179 PetscScalar one = 1.0; 180 int mX,mY,nX,nY,ierr; 181 182 PetscFunctionBegin; 183 PetscValidHeaderSpecific(X,MAT_COOKIE); 184 PetscValidHeaderSpecific(Y,MAT_COOKIE); 185 PetscValidScalarPointer(a); 186 187 ierr = MatGetSize(X,&mX,&nX);CHKERRQ(ierr); 188 ierr = MatGetSize(X,&mY,&nY);CHKERRQ(ierr); 189 if (mX != mY || nX != nY) SETERRQ4(PETSC_ERR_ARG_SIZ,"Non conforming matrices: %d %d first %d %d second",mX,mY,nX,nY); 190 191 ierr = MatScale(a,Y);CHKERRQ(ierr); 192 ierr = MatAXPY(&one,X,Y);CHKERRQ(ierr); 193 PetscFunctionReturn(0); 194 } 195 196 #undef __FUNCT__ 197 #define __FUNCT__ "MatComputeExplicitOperator" 198 /*@ 199 MatComputeExplicitOperator - Computes the explicit matrix 200 201 Collective on Mat 202 203 Input Parameter: 204 . inmat - the matrix 205 206 Output Parameter: 207 . mat - the explict preconditioned operator 208 209 Notes: 210 This computation is done by applying the operators to columns of the 211 identity matrix. 212 213 Currently, this routine uses a dense matrix format when 1 processor 214 is used and a sparse format otherwise. This routine is costly in general, 215 and is recommended for use only with relatively small systems. 216 217 Level: advanced 218 219 .keywords: Mat, compute, explicit, operator 220 221 @*/ 222 int MatComputeExplicitOperator(Mat inmat,Mat *mat) 223 { 224 Vec in,out; 225 int ierr,i,M,m,size,*rows,start,end; 226 MPI_Comm comm; 227 PetscScalar *array,zero = 0.0,one = 1.0; 228 229 PetscFunctionBegin; 230 PetscValidHeaderSpecific(inmat,MAT_COOKIE); 231 PetscValidPointer(mat); 232 233 comm = inmat->comm; 234 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 235 236 ierr = MatGetLocalSize(inmat,&m,0);CHKERRQ(ierr); 237 ierr = MatGetSize(inmat,&M,0);CHKERRQ(ierr); 238 ierr = VecCreateMPI(comm,m,M,&in);CHKERRQ(ierr); 239 ierr = VecDuplicate(in,&out);CHKERRQ(ierr); 240 ierr = VecGetOwnershipRange(in,&start,&end);CHKERRQ(ierr); 241 ierr = PetscMalloc((m+1)*sizeof(int),&rows);CHKERRQ(ierr); 242 for (i=0; i<m; i++) {rows[i] = start + i;} 243 244 if (size == 1) { 245 ierr = MatCreateSeqDense(comm,M,M,PETSC_NULL,mat);CHKERRQ(ierr); 246 } else { 247 ierr = MatCreateMPIAIJ(comm,m,m,M,M,0,0,0,0,mat);CHKERRQ(ierr); 248 } 249 250 for (i=0; i<M; i++) { 251 252 ierr = VecSet(&zero,in);CHKERRQ(ierr); 253 ierr = VecSetValues(in,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr); 254 ierr = VecAssemblyBegin(in);CHKERRQ(ierr); 255 ierr = VecAssemblyEnd(in);CHKERRQ(ierr); 256 257 ierr = MatMult(inmat,in,out);CHKERRQ(ierr); 258 259 ierr = VecGetArray(out,&array);CHKERRQ(ierr); 260 ierr = MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 261 ierr = VecRestoreArray(out,&array);CHKERRQ(ierr); 262 263 } 264 ierr = PetscFree(rows);CHKERRQ(ierr); 265 ierr = VecDestroy(out);CHKERRQ(ierr); 266 ierr = VecDestroy(in);CHKERRQ(ierr); 267 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 268 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 272