1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: axpy.c,v 1.28 1997/09/11 20:40:16 bsmith Exp bsmith $"; 3 #endif 4 5 #include "src/mat/matimpl.h" /*I "mat.h" I*/ 6 7 #undef __FUNC__ 8 #define __FUNC__ "MatAXPY" 9 /*@ 10 MatAXPY - Computes Y = a*X + Y. 11 12 Input Parameters: 13 . X,Y - the matrices 14 . a - the scalar multiplier 15 16 Contributed by: Matthew Knepley 17 18 .keywords: matrix, add 19 20 @*/ 21 int MatAXPY(Scalar *a,Mat X,Mat Y) 22 { 23 int m1,m2,n1,n2,i,*row,start,end,j,ncols,ierr; 24 Scalar *val,*vals; 25 26 PetscFunctionBegin; 27 PetscValidHeaderSpecific(X,MAT_COOKIE); 28 PetscValidHeaderSpecific(Y,MAT_COOKIE); 29 PetscValidScalarPointer(a); 30 31 MatGetSize(X,&m1,&n1); MatGetSize(Y,&m2,&n2); 32 if (m1 != m2 || n1 != n2) SETERRQ(1,0,"Non conforming matrix add"); 33 34 if (X->ops.axpy) { 35 ierr = (*X->ops.axpy)(a,X,Y); CHKERRQ(ierr); 36 } else { 37 ierr = MatGetOwnershipRange(X,&start,&end); CHKERRQ(ierr); 38 if (*a == 1.0) { 39 for (i = start; i < end; i++) { 40 ierr = MatGetRow(X, i, &ncols, &row, &vals); CHKERRQ(ierr); 41 ierr = MatSetValues(Y, 1, &i, ncols, row, vals, ADD_VALUES); CHKERRQ(ierr); 42 ierr = MatRestoreRow(X, i, &ncols, &row, &vals); CHKERRQ(ierr); 43 } 44 } else { 45 vals = (Scalar *) PetscMalloc( (n1+1)*sizeof(Scalar) ); CHKPTRQ(vals); 46 for ( i=start; i<end; i++ ) { 47 ierr = MatGetRow(X,i,&ncols,&row,&val); CHKERRQ(ierr); 48 for ( j=0; j<ncols; j++ ) { 49 vals[j] = (*a)*val[j]; 50 } 51 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES); CHKERRQ(ierr); 52 ierr = MatRestoreRow(X,i,&ncols,&row,&val); CHKERRQ(ierr); 53 } 54 PetscFree(vals); 55 } 56 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 57 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 58 } 59 PetscFunctionReturn(0); 60 } 61 62 #undef __FUNC__ 63 #define __FUNC__ "MatShift" 64 /*@ 65 MatShift - Computes Y = Y + a I, where a is a scalar and I is the identity 66 matrix. 67 68 Input Parameters: 69 . Y - the matrices 70 . a - the scalar 71 72 .keywords: matrix, add, shift 73 74 .seealso: MatDiagonalShift() 75 @*/ 76 int MatShift(Scalar *a,Mat Y) 77 { 78 int i,start,end,ierr; 79 80 PetscFunctionBegin; 81 PetscValidHeaderSpecific(Y,MAT_COOKIE); 82 PetscValidScalarPointer(a); 83 if (Y->ops.shift) { 84 ierr = (*Y->ops.shift)(a,Y); CHKERRQ(ierr); 85 } 86 else { 87 ierr = MatGetOwnershipRange(Y,&start,&end); CHKERRQ(ierr); 88 for ( i=start; i<end; i++ ) { 89 ierr = MatSetValues(Y,1,&i,1,&i,a,ADD_VALUES); CHKERRQ(ierr); 90 } 91 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 92 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 93 } 94 PetscFunctionReturn(0); 95 } 96 97 #undef __FUNC__ 98 #define __FUNC__ "MatDiagonalShift" 99 /*@ 100 MatDiagonalShift - Computes Y = Y + D, where D is a diagonal matrix 101 that is represented as a vector. 102 103 Input Parameters: 104 . Y - the input matrix 105 . D - the diagonal matrix, represented as a vector 106 107 Input Parameters: 108 . Y - the shifted ouput matrix 109 110 .keywords: matrix, add, shift, diagonal 111 112 .seealso: MatShift() 113 @*/ 114 int MatDiagonalShift(Mat Y,Vec D) 115 { 116 int i,start,end,ierr; 117 118 PetscFunctionBegin; 119 PetscValidHeaderSpecific(Y,MAT_COOKIE); 120 PetscValidHeaderSpecific(D,VEC_COOKIE); 121 if (Y->ops.shift) { 122 ierr = (*Y->ops.diagonalshift)(D,Y); CHKERRQ(ierr); 123 } 124 else { 125 int vstart,vend; 126 Scalar *v; 127 ierr = VecGetOwnershipRange(D,&vstart,&vend); CHKERRQ(ierr); 128 ierr = MatGetOwnershipRange(Y,&start,&end); CHKERRQ(ierr); 129 if (vstart != start || vend != end) { 130 SETERRQ(1,0,"Vector ownership range not compatible with matrix"); 131 } 132 133 ierr = VecGetArray(D,&v); CHKERRQ(ierr); 134 for ( i=start; i<end; i++ ) { 135 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,ADD_VALUES); CHKERRQ(ierr); 136 } 137 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 138 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 139 } 140 PetscFunctionReturn(0); 141 } 142 143 #undef __FUNC__ 144 #define __FUNC__ "MatAYPX" 145 /*@ 146 MatAYPX - Computes Y = X + a*Y. 147 148 Input Parameters: 149 . X,Y - the matrices 150 . a - the scalar multiplier 151 152 Contributed by: Matthew Knepley 153 154 .keywords: matrix, add 155 156 @*/ 157 int MatAYPX(Scalar *a,Mat X,Mat Y) 158 { 159 Scalar one = 1.0; 160 int mX, mY,nX, nY,ierr; 161 162 PetscFunctionBegin; 163 PetscValidHeaderSpecific(X, MAT_COOKIE); 164 PetscValidHeaderSpecific(Y, MAT_COOKIE); 165 PetscValidScalarPointer(a); 166 167 MatGetSize(X, &mX, &nX); 168 MatGetSize(X, &mY, &nY); 169 if (mX != mY || nX != nY) SETERRQ(1,0,"Non conforming matrices"); 170 171 ierr = MatScale(a, Y); CHKERRQ(ierr); 172 ierr = MatAXPY(&one, X, Y); CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175