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