1 /*$Id: axpy.c,v 1.54 2001/08/06 21:16:10 bsmith Exp bsmith $*/ 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 + a - the scalar multiplier 14 . X - the first matrix 15 . Y - the second matrix 16 - str - either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN 17 18 Contributed by: Matthew Knepley 19 20 Notes: 21 Will only be efficient if one has the SAME_NONZERO_PATTERN or SUBSET_NONZERO_PATTERN 22 23 Level: intermediate 24 25 .keywords: matrix, add 26 27 .seealso: MatAYPX() 28 @*/ 29 int MatAXPY(const PetscScalar *a,Mat X,Mat Y,MatStructure str) 30 { 31 int m1,m2,n1,n2,ierr; 32 33 PetscFunctionBegin; 34 PetscValidScalarPointer(a,1); 35 PetscValidHeaderSpecific(X,MAT_COOKIE,2); 36 PetscValidHeaderSpecific(Y,MAT_COOKIE,3); 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,str);CHKERRQ(ierr); 44 } else { 45 ierr = MatAXPY_Basic(a,X,Y,str);CHKERRQ(ierr); 46 } 47 PetscFunctionReturn(0); 48 } 49 50 51 #undef __FUNCT__ 52 #define __FUNCT__ "MatAXPY_Basic" 53 int MatAXPY_Basic(const PetscScalar *a,Mat X,Mat Y,MatStructure str) 54 { 55 int i,*row,start,end,j,ncols,ierr,m,n; 56 PetscScalar *val,*vals; 57 58 PetscFunctionBegin; 59 ierr = MatGetSize(X,&m,&n);CHKERRQ(ierr); 60 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 61 if (*a == 1.0) { 62 for (i = start; i < end; i++) { 63 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 64 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 65 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 66 } 67 } else { 68 ierr = PetscMalloc((n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 69 for (i=start; i<end; i++) { 70 ierr = MatGetRow(X,i,&ncols,&row,&val);CHKERRQ(ierr); 71 for (j=0; j<ncols; j++) { 72 vals[j] = (*a)*val[j]; 73 } 74 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 75 ierr = MatRestoreRow(X,i,&ncols,&row,&val);CHKERRQ(ierr); 76 } 77 ierr = PetscFree(vals);CHKERRQ(ierr); 78 } 79 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 80 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "MatShift" 86 /*@ 87 MatShift - Computes Y = Y + a I, where a is a PetscScalar and I is the identity matrix. 88 89 Collective on Mat 90 91 Input Parameters: 92 + Y - the matrices 93 - a - the PetscScalar 94 95 Level: intermediate 96 97 .keywords: matrix, add, shift 98 99 .seealso: MatDiagonalSet() 100 @*/ 101 int MatShift(const PetscScalar *a,Mat Y) 102 { 103 int i,start,end,ierr; 104 105 PetscFunctionBegin; 106 PetscValidScalarPointer(a,1); 107 PetscValidHeaderSpecific(Y,MAT_COOKIE,2); 108 if (Y->ops->shift) { 109 ierr = (*Y->ops->shift)(a,Y);CHKERRQ(ierr); 110 } else { 111 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 112 for (i=start; i<end; i++) { 113 ierr = MatSetValues(Y,1,&i,1,&i,a,ADD_VALUES);CHKERRQ(ierr); 114 } 115 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 116 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 117 } 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNCT__ 122 #define __FUNCT__ "MatDiagonalSet" 123 /*@ 124 MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix 125 that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is 126 INSERT_VALUES. 127 128 Input Parameters: 129 + Y - the input matrix 130 . D - the diagonal matrix, represented as a vector 131 - i - INSERT_VALUES or ADD_VALUES 132 133 Collective on Mat and Vec 134 135 Level: intermediate 136 137 .keywords: matrix, add, shift, diagonal 138 139 .seealso: MatShift() 140 @*/ 141 int MatDiagonalSet(Mat Y,Vec D,InsertMode is) 142 { 143 int i,start,end,ierr; 144 145 PetscFunctionBegin; 146 PetscValidHeaderSpecific(Y,MAT_COOKIE,1); 147 PetscValidHeaderSpecific(D,VEC_COOKIE,2); 148 if (Y->ops->diagonalset) { 149 ierr = (*Y->ops->diagonalset)(Y,D,is);CHKERRQ(ierr); 150 } else { 151 int vstart,vend; 152 PetscScalar *v; 153 ierr = VecGetOwnershipRange(D,&vstart,&vend);CHKERRQ(ierr); 154 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 155 if (vstart != start || vend != end) { 156 SETERRQ4(PETSC_ERR_ARG_SIZ,"Vector ownership range not compatible with matrix: %d %d vec %d %d mat",vstart,vend,start,end); 157 } 158 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 159 for (i=start; i<end; i++) { 160 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,is);CHKERRQ(ierr); 161 } 162 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 163 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 164 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 165 } 166 PetscFunctionReturn(0); 167 } 168 169 #undef __FUNCT__ 170 #define __FUNCT__ "MatAYPX" 171 /*@ 172 MatAYPX - Computes Y = X + a*Y. 173 174 Collective on Mat 175 176 Input Parameters: 177 + X,Y - the matrices 178 - a - the PetscScalar multiplier 179 180 Contributed by: Matthew Knepley 181 182 Notes: 183 This routine currently uses the MatAXPY() implementation. 184 185 This is slow, if you need it fast send email to petsc-maint@mcs.anl.gov 186 187 Level: intermediate 188 189 .keywords: matrix, add 190 191 .seealso: MatAXPY() 192 @*/ 193 int MatAYPX(const PetscScalar *a,Mat X,Mat Y) 194 { 195 PetscScalar one = 1.0; 196 int mX,mY,nX,nY,ierr; 197 198 PetscFunctionBegin; 199 PetscValidScalarPointer(a,1); 200 PetscValidHeaderSpecific(X,MAT_COOKIE,2); 201 PetscValidHeaderSpecific(Y,MAT_COOKIE,3); 202 203 ierr = MatGetSize(X,&mX,&nX);CHKERRQ(ierr); 204 ierr = MatGetSize(X,&mY,&nY);CHKERRQ(ierr); 205 if (mX != mY || nX != nY) SETERRQ4(PETSC_ERR_ARG_SIZ,"Non conforming matrices: %d %d first %d %d second",mX,mY,nX,nY); 206 207 ierr = MatScale(a,Y);CHKERRQ(ierr); 208 ierr = MatAXPY(&one,X,Y,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 209 PetscFunctionReturn(0); 210 } 211 212 #undef __FUNCT__ 213 #define __FUNCT__ "MatComputeExplicitOperator" 214 /*@ 215 MatComputeExplicitOperator - Computes the explicit matrix 216 217 Collective on Mat 218 219 Input Parameter: 220 . inmat - the matrix 221 222 Output Parameter: 223 . mat - the explict preconditioned operator 224 225 Notes: 226 This computation is done by applying the operators to columns of the 227 identity matrix. 228 229 Currently, this routine uses a dense matrix format when 1 processor 230 is used and a sparse format otherwise. This routine is costly in general, 231 and is recommended for use only with relatively small systems. 232 233 Level: advanced 234 235 .keywords: Mat, compute, explicit, operator 236 237 @*/ 238 int MatComputeExplicitOperator(Mat inmat,Mat *mat) 239 { 240 Vec in,out; 241 int ierr,i,M,m,size,*rows,start,end; 242 MPI_Comm comm; 243 PetscScalar *array,zero = 0.0,one = 1.0; 244 245 PetscFunctionBegin; 246 PetscValidHeaderSpecific(inmat,MAT_COOKIE,1); 247 PetscValidPointer(mat,2); 248 249 comm = inmat->comm; 250 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 251 252 ierr = MatGetLocalSize(inmat,&m,0);CHKERRQ(ierr); 253 ierr = MatGetSize(inmat,&M,0);CHKERRQ(ierr); 254 ierr = VecCreateMPI(comm,m,M,&in);CHKERRQ(ierr); 255 ierr = VecDuplicate(in,&out);CHKERRQ(ierr); 256 ierr = VecGetOwnershipRange(in,&start,&end);CHKERRQ(ierr); 257 ierr = PetscMalloc((m+1)*sizeof(int),&rows);CHKERRQ(ierr); 258 for (i=0; i<m; i++) {rows[i] = start + i;} 259 260 if (size == 1) { 261 ierr = MatCreateSeqDense(comm,M,M,PETSC_NULL,mat);CHKERRQ(ierr); 262 } else { 263 ierr = MatCreateMPIAIJ(comm,m,m,M,M,0,0,0,0,mat);CHKERRQ(ierr); 264 } 265 266 for (i=0; i<M; i++) { 267 268 ierr = VecSet(&zero,in);CHKERRQ(ierr); 269 ierr = VecSetValues(in,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr); 270 ierr = VecAssemblyBegin(in);CHKERRQ(ierr); 271 ierr = VecAssemblyEnd(in);CHKERRQ(ierr); 272 273 ierr = MatMult(inmat,in,out);CHKERRQ(ierr); 274 275 ierr = VecGetArray(out,&array);CHKERRQ(ierr); 276 ierr = MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 277 ierr = VecRestoreArray(out,&array);CHKERRQ(ierr); 278 279 } 280 ierr = PetscFree(rows);CHKERRQ(ierr); 281 ierr = VecDestroy(out);CHKERRQ(ierr); 282 ierr = VecDestroy(in);CHKERRQ(ierr); 283 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 284 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 285 PetscFunctionReturn(0); 286 } 287 288 /* Get the map xtoy which is used by MatAXPY() in the case of SUBSET_NONZERO_PATTERN */ 289 #undef __FUNCT__ 290 #define __FUNCT__ "MatAXPYGetxtoy_Private" 291 int MatAXPYGetxtoy_Private(int m,int *xi,int *xj,int *xgarray, int *yi,int *yj,int *ygarray, int **xtoy) 292 { 293 int ierr,row,i,nz,xcol,ycol,jx,jy,*x2y; 294 295 PetscFunctionBegin; 296 ierr = PetscMalloc(xi[m]*sizeof(int),&x2y);CHKERRQ(ierr); 297 i = 0; 298 for (row=0; row<m; row++){ 299 nz = xi[1] - xi[0]; 300 jy = 0; 301 for (jx=0; jx<nz; jx++,jy++){ 302 if (xgarray && ygarray){ 303 xcol = xgarray[xj[*xi + jx]]; 304 ycol = ygarray[yj[*yi + jy]]; 305 } else { 306 xcol = xj[*xi + jx]; 307 ycol = yj[*yi + jy]; /* col index for y */ 308 } 309 while ( ycol < xcol ) { 310 jy++; 311 if (ygarray){ 312 ycol = ygarray[yj[*yi + jy]]; 313 } else { 314 ycol = yj[*yi + jy]; 315 } 316 } 317 if (xcol != ycol) SETERRQ2(PETSC_ERR_ARG_WRONG,"X matrix entry (%d,%d) is not in Y matrix",row,ycol); 318 x2y[i++] = *yi + jy; 319 } 320 xi++; yi++; 321 } 322 *xtoy = x2y; 323 PetscFunctionReturn(0); 324 } 325