1 /*$Id: axpy.c,v 1.48 2000/09/28 21:12:12 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 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 vals = (Scalar*)PetscMalloc((n1+1)*sizeof(Scalar));CHKPTRQ(vals); 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 __FUNC__ 71 #define __FUNC__ /*<a name=""></a>*/"MatShift" 72 /*@ 73 MatShift - Computes Y = Y + a I, where a is a scalar and I is the identity matrix. 74 75 Collective on Mat 76 77 Input Parameters: 78 + Y - the matrices 79 - a - the scalar 80 81 Level: intermediate 82 83 .keywords: matrix, add, shift 84 85 .seealso: MatDiagonalSet() 86 @*/ 87 int MatShift(Scalar *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 __FUNC__ 108 #define __FUNC__ /*<a name="MatDiagonalSet"></a>*/"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 Scalar *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 __FUNC__ 156 #define __FUNC__ /*<a name=""></a>*/"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 scalar 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(Scalar *a,Mat X,Mat Y) 178 { 179 Scalar 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