1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: axpy.c,v 1.32 1998/04/13 17:43:46 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 } 90 else { 91 ierr = MatGetOwnershipRange(Y,&start,&end); CHKERRQ(ierr); 92 for ( i=start; i<end; i++ ) { 93 ierr = MatSetValues(Y,1,&i,1,&i,a,ADD_VALUES); CHKERRQ(ierr); 94 } 95 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 96 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 97 } 98 PetscFunctionReturn(0); 99 } 100 101 #undef __FUNC__ 102 #define __FUNC__ "MatDiagonalShift" 103 /*@ 104 MatDiagonalShift - Computes Y = Y + D, where D is a diagonal matrix 105 that is represented as a vector. 106 107 Input Parameters: 108 + Y - the input matrix 109 - D - the diagonal matrix, represented as a vector 110 111 Input Parameters: 112 . Y - the shifted ouput matrix 113 114 Collective on Mat and Vec 115 116 .keywords: matrix, add, shift, diagonal 117 118 .seealso: MatShift() 119 @*/ 120 int MatDiagonalShift(Mat Y,Vec D) 121 { 122 int i,start,end,ierr; 123 124 PetscFunctionBegin; 125 PetscValidHeaderSpecific(Y,MAT_COOKIE); 126 PetscValidHeaderSpecific(D,VEC_COOKIE); 127 if (Y->ops->shift) { 128 ierr = (*Y->ops->diagonalshift)(D,Y); CHKERRQ(ierr); 129 } 130 else { 131 int vstart,vend; 132 Scalar *v; 133 ierr = VecGetOwnershipRange(D,&vstart,&vend); CHKERRQ(ierr); 134 ierr = MatGetOwnershipRange(Y,&start,&end); CHKERRQ(ierr); 135 if (vstart != start || vend != end) { 136 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Vector ownership range not compatible with matrix"); 137 } 138 139 ierr = VecGetArray(D,&v); CHKERRQ(ierr); 140 for ( i=start; i<end; i++ ) { 141 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,ADD_VALUES); CHKERRQ(ierr); 142 } 143 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 144 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 145 } 146 PetscFunctionReturn(0); 147 } 148 149 #undef __FUNC__ 150 #define __FUNC__ "MatAYPX" 151 /*@ 152 MatAYPX - Computes Y = X + a*Y. 153 154 Collective on Mat 155 156 Input Parameters: 157 + X,Y - the matrices 158 - a - the scalar multiplier 159 160 Contributed by: Matthew Knepley 161 162 .keywords: matrix, add 163 164 @*/ 165 int MatAYPX(Scalar *a,Mat X,Mat Y) 166 { 167 Scalar one = 1.0; 168 int mX, mY,nX, nY,ierr; 169 170 PetscFunctionBegin; 171 PetscValidHeaderSpecific(X, MAT_COOKIE); 172 PetscValidHeaderSpecific(Y, MAT_COOKIE); 173 PetscValidScalarPointer(a); 174 175 MatGetSize(X, &mX, &nX); 176 MatGetSize(X, &mY, &nY); 177 if (mX != mY || nX != nY) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Non conforming matrices"); 178 179 ierr = MatScale(a, Y); CHKERRQ(ierr); 180 ierr = MatAXPY(&one, X, Y); CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183