1 /*$Id: axpy.c,v 1.47 2000/05/28 15:21:54 bsmith Exp bsmith $*/ 2 3 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 4 5 #undef __FUNC__ 6 #define __FUNC__ /*<a name=""></a>*/"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 scalar 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(Scalar *a,Mat X,Mat Y) 29 { 30 int m1,m2,n1,n2,i,*row,start,end,j,ncols,ierr; 31 Scalar *val,*vals; 32 33 PetscFunctionBegin; 34 PetscValidHeaderSpecific(X,MAT_COOKIE); 35 PetscValidHeaderSpecific(Y,MAT_COOKIE); 36 PetscValidScalarPointer(a); 37 38 MatGetSize(X,&m1,&n1); MatGetSize(Y,&m2,&n2); 39 if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_ERR_ARG_SIZ,"Non conforming matrix add: %d %d %d %d",m1,m2,n1,n2); 40 41 if (X->ops->axpy) { 42 ierr = (*X->ops->axpy)(a,X,Y);CHKERRQ(ierr); 43 } else { 44 ierr = MatGetOwnershipRange(X,&start,&end);CHKERRQ(ierr); 45 if (*a == 1.0) { 46 for (i = start; i < end; i++) { 47 ierr = MatGetRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 48 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 49 ierr = MatRestoreRow(X,i,&ncols,&row,&vals);CHKERRQ(ierr); 50 } 51 } else { 52 vals = (Scalar*)PetscMalloc((n1+1)*sizeof(Scalar));CHKPTRQ(vals); 53 for (i=start; i<end; i++) { 54 ierr = MatGetRow(X,i,&ncols,&row,&val);CHKERRQ(ierr); 55 for (j=0; j<ncols; j++) { 56 vals[j] = (*a)*val[j]; 57 } 58 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES);CHKERRQ(ierr); 59 ierr = MatRestoreRow(X,i,&ncols,&row,&val);CHKERRQ(ierr); 60 } 61 ierr = PetscFree(vals);CHKERRQ(ierr); 62 } 63 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 64 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 65 } 66 PetscFunctionReturn(0); 67 } 68 69 #undef __FUNC__ 70 #define __FUNC__ /*<a name=""></a>*/"MatShift" 71 /*@ 72 MatShift - Computes Y = Y + a I, where a is a scalar and I is the identity matrix. 73 74 Collective on Mat 75 76 Input Parameters: 77 + Y - the matrices 78 - a - the scalar 79 80 Level: intermediate 81 82 .keywords: matrix, add, shift 83 84 .seealso: MatDiagonalSet() 85 @*/ 86 int MatShift(Scalar *a,Mat Y) 87 { 88 int i,start,end,ierr; 89 90 PetscFunctionBegin; 91 PetscValidHeaderSpecific(Y,MAT_COOKIE); 92 PetscValidScalarPointer(a); 93 if (Y->ops->shift) { 94 ierr = (*Y->ops->shift)(a,Y);CHKERRQ(ierr); 95 } else { 96 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 97 for (i=start; i<end; i++) { 98 ierr = MatSetValues(Y,1,&i,1,&i,a,ADD_VALUES);CHKERRQ(ierr); 99 } 100 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 101 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 102 } 103 PetscFunctionReturn(0); 104 } 105 106 #undef __FUNC__ 107 #define __FUNC__ /*<a name="MatDiagonalSet"></a>*/"MatDiagonalSet" 108 /*@ 109 MatDiagonalSet - Computes Y = Y + D, where D is a diagonal matrix 110 that is represented as a vector. Or Y[i,i] = D[i] if InsertMode is 111 INSERT_VALUES. 112 113 Input Parameters: 114 + Y - the input matrix 115 . D - the diagonal matrix, represented as a vector 116 - i - INSERT_VALUES or ADD_VALUES 117 118 Collective on Mat and Vec 119 120 Level: intermediate 121 122 .keywords: matrix, add, shift, diagonal 123 124 .seealso: MatShift() 125 @*/ 126 int MatDiagonalSet(Mat Y,Vec D,InsertMode is) 127 { 128 int i,start,end,ierr; 129 130 PetscFunctionBegin; 131 PetscValidHeaderSpecific(Y,MAT_COOKIE); 132 PetscValidHeaderSpecific(D,VEC_COOKIE); 133 if (Y->ops->diagonalset) { 134 ierr = (*Y->ops->diagonalset)(Y,D,is);CHKERRQ(ierr); 135 } else { 136 int vstart,vend; 137 Scalar *v; 138 ierr = VecGetOwnershipRange(D,&vstart,&vend);CHKERRQ(ierr); 139 ierr = MatGetOwnershipRange(Y,&start,&end);CHKERRQ(ierr); 140 if (vstart != start || vend != end) { 141 SETERRQ4(PETSC_ERR_ARG_SIZ,"Vector ownership range not compatible with matrix: %d %d vec %d %d mat",vstart,vend,start,end); 142 } 143 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 144 for (i=start; i<end; i++) { 145 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,is);CHKERRQ(ierr); 146 } 147 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 148 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 149 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 150 } 151 PetscFunctionReturn(0); 152 } 153 154 #undef __FUNC__ 155 #define __FUNC__ /*<a name=""></a>*/"MatAYPX" 156 /*@ 157 MatAYPX - Computes Y = X + a*Y. 158 159 Collective on Mat 160 161 Input Parameters: 162 + X,Y - the matrices 163 - a - the scalar multiplier 164 165 Contributed by: Matthew Knepley 166 167 Notes: 168 This routine currently uses the MatAXPY() implementation. 169 170 Level: intermediate 171 172 .keywords: matrix, add 173 174 .seealso: MatAXPY() 175 @*/ 176 int MatAYPX(Scalar *a,Mat X,Mat Y) 177 { 178 Scalar one = 1.0; 179 int mX,mY,nX,nY,ierr; 180 181 PetscFunctionBegin; 182 PetscValidHeaderSpecific(X,MAT_COOKIE); 183 PetscValidHeaderSpecific(Y,MAT_COOKIE); 184 PetscValidScalarPointer(a); 185 186 ierr = MatGetSize(X,&mX,&nX);CHKERRQ(ierr); 187 ierr = MatGetSize(X,&mY,&nY);CHKERRQ(ierr); 188 if (mX != mY || nX != nY) SETERRQ4(PETSC_ERR_ARG_SIZ,"Non conforming matrices: %d %d first %d %d second",mX,mY,nX,nY); 189 190 ierr = MatScale(a,Y);CHKERRQ(ierr); 191 ierr = MatAXPY(&one,X,Y);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194