1 /*$Id: axpy.c,v 1.43 2000/04/09 04:36:57 bsmith Exp bsmith $*/ 2 3 #include "src/mat/matimpl.h" /*I "mat.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,0,"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: MatDiagonalShift() 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=""></a>*/"MatDiagonalShift" 108 /*@ 109 MatDiagonalShift - Computes Y = Y + D, where D is a diagonal matrix 110 that is represented as a vector. 111 112 Input Parameters: 113 + Y - the input matrix 114 - D - the diagonal matrix, represented as a vector 115 116 Input Parameters: 117 . Y - the shifted ouput matrix 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 MatDiagonalShift(Mat Y,Vec D) 128 { 129 int i,start,end,ierr; 130 131 PetscFunctionBegin; 132 PetscValidHeaderSpecific(Y,MAT_COOKIE); 133 PetscValidHeaderSpecific(D,VEC_COOKIE); 134 if (Y->ops->shift) { 135 ierr = (*Y->ops->diagonalshift)(D,Y);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,0,"Vector ownership range not compatible with matrix: %d %d vec %d %d mat",vstart,vend,start,end); 143 } 144 145 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 146 for (i=start; i<end; i++) { 147 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,ADD_VALUES);CHKERRQ(ierr); 148 } 149 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 150 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 151 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 152 } 153 PetscFunctionReturn(0); 154 } 155 156 #undef __FUNC__ 157 #define __FUNC__ /*<a name=""></a>*/"MatAYPX" 158 /*@ 159 MatAYPX - Computes Y = X + a*Y. 160 161 Collective on Mat 162 163 Input Parameters: 164 + X,Y - the matrices 165 - a - the scalar multiplier 166 167 Contributed by: Matthew Knepley 168 169 Notes: 170 This routine currently uses the MatAXPY() implementation. 171 172 Level: intermediate 173 174 .keywords: matrix, add 175 176 .seealso: MatAXPY() 177 @*/ 178 int MatAYPX(Scalar *a,Mat X,Mat Y) 179 { 180 Scalar one = 1.0; 181 int mX,mY,nX,nY,ierr; 182 183 PetscFunctionBegin; 184 PetscValidHeaderSpecific(X,MAT_COOKIE); 185 PetscValidHeaderSpecific(Y,MAT_COOKIE); 186 PetscValidScalarPointer(a); 187 188 ierr = MatGetSize(X,&mX,&nX);CHKERRQ(ierr); 189 ierr = MatGetSize(X,&mY,&nY);CHKERRQ(ierr); 190 if (mX != mY || nX != nY) SETERRQ4(PETSC_ERR_ARG_SIZ,0,"Non conforming matrices: %d %d first %d %d second",mX,mY,nX,nY); 191 192 ierr = MatScale(a,Y);CHKERRQ(ierr); 193 ierr = MatAXPY(&one,X,Y);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196