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