1 #ifndef lint 2 static char vcid[] = "$Id: axpy.c,v 1.24 1997/01/06 20:26:03 balay Exp balay $"; 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 .keywords: matrix, add 17 @*/ 18 int MatAXPY(Scalar *a,Mat X,Mat Y) 19 { 20 int m1,m2,n1,n2,i,*row,start,end,j,ncols,ierr; 21 Scalar *val,*vals; 22 23 PetscValidHeaderSpecific(X,MAT_COOKIE); PetscValidHeaderSpecific(Y,MAT_COOKIE); 24 PetscValidScalarPointer(a); 25 26 MatGetSize(X,&m1,&n1); MatGetSize(X,&m2,&n2); 27 if (m1 != m2 || n1 != n2) SETERRQ(1,0,"Non conforming matrix add"); 28 29 if (X->ops.axpy) { 30 ierr = (*X->ops.axpy)(a,X,Y); CHKERRQ(ierr); 31 } 32 else { 33 vals = (Scalar *) PetscMalloc( (n1+1)*sizeof(Scalar) ); CHKPTRQ(vals); 34 ierr = MatGetOwnershipRange(X,&start,&end); CHKERRQ(ierr); 35 for ( i=start; i<end; i++ ) { 36 ierr = MatGetRow(X,i,&ncols,&row,&val); CHKERRQ(ierr); 37 for ( j=0; j<ncols; j++ ) { 38 vals[j] = (*a)*val[j]; 39 } 40 ierr = MatSetValues(Y,1,&i,ncols,row,vals,ADD_VALUES); CHKERRQ(ierr); 41 ierr = MatRestoreRow(X,i,&ncols,&row,&val); CHKERRQ(ierr); 42 } 43 PetscFree(vals); 44 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 45 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 46 } 47 return 0; 48 } 49 50 #undef __FUNC__ 51 #define __FUNC__ "MatShift" 52 /*@ 53 MatShift - Computes Y = Y + a I, where a is a scalar and I is the identity 54 matrix. 55 56 Input Parameters: 57 . Y - the matrices 58 . a - the scalar 59 60 .keywords: matrix, add, shift 61 62 .seealso: MatDiagonalShift() 63 @*/ 64 int MatShift(Scalar *a,Mat Y) 65 { 66 int i,start,end,ierr; 67 68 PetscValidHeaderSpecific(Y,MAT_COOKIE); 69 PetscValidScalarPointer(a); 70 if (Y->ops.shift) { 71 ierr = (*Y->ops.shift)(a,Y); CHKERRQ(ierr); 72 } 73 else { 74 MatGetOwnershipRange(Y,&start,&end); 75 for ( i=start; i<end; i++ ) { 76 ierr = MatSetValues(Y,1,&i,1,&i,a,ADD_VALUES); CHKERRQ(ierr); 77 } 78 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 79 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 80 } 81 return 0; 82 } 83 84 #undef __FUNC__ 85 #define __FUNC__ "MatDiagonalShift" 86 /*@ 87 MatDiagonalShift - Computes Y = Y + D, where D is a diagonal matrix 88 that is represented as a vector. 89 90 Input Parameters: 91 . Y - the input matrix 92 . D - the diagonal matrix, represented as a vector 93 94 Input Parameters: 95 . Y - the shifted ouput matrix 96 97 .keywords: matrix, add, shift, diagonal 98 99 .seealso: MatShift() 100 @*/ 101 int MatDiagonalShift(Mat Y,Vec D) 102 { 103 int i,start,end,ierr; 104 105 PetscValidHeaderSpecific(Y,MAT_COOKIE); 106 PetscValidHeaderSpecific(D,VEC_COOKIE); 107 if (Y->ops.shift) { 108 ierr = (*Y->ops.diagonalshift)(D,Y); CHKERRQ(ierr); 109 } 110 else { 111 int vstart,vend; 112 Scalar *v; 113 ierr = VecGetOwnershipRange(D,&vstart,&vend); CHKERRQ(ierr); 114 ierr = MatGetOwnershipRange(Y,&start,&end); CHKERRQ(ierr); 115 if (vstart != start || vend != end) 116 SETERRQ(1,0,"Vector shift not compatible with matrix"); 117 118 ierr = VecGetArray(D,&v); CHKERRQ(ierr); 119 for ( i=start; i<end; i++ ) { 120 ierr = MatSetValues(Y,1,&i,1,&i,v+i-start,ADD_VALUES); CHKERRQ(ierr); 121 } 122 ierr = MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 123 ierr = MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 124 } 125 return 0; 126 } 127