1 subroutine SparseAp(iper, ilwork, iBC, col, row, lhsK, p) 2 3C============================================================================ 4C 5C "SparseAp": This routine performs the product of a sparsley stored matrix 6C lhsK(nflow*nflow, nnz_tot) and a vector p(nshg, nflow). The 7C results of the product is returned in p(nshg, nflow). 8C 9C Nahid Razmara, Spring 2000. (Sparse Matrix) 10C============================================================================ 11 12c 13 include "common.h" 14c 15c 16c.... Data declaration 17c 18 integer col(nshg+1), row(nnz*nshg), 19 & iper(nshg), ilwork(nlwork) 20 real*8 lhsK(nflow*nflow,nnz_tot) 21 real*8 p(nshg,nflow), q(nshg,nflow) 22c 23 real*8 tmp1, tmp2, tmp3, tmp4, tmp5 24c 25c.... communicate:: copy the master's portion of uBrg to each slave 26c 27 if (numpe > 1) then 28 call commu (p, ilwork, nflow , 'out') 29 endif 30c 31c.... local periodic boundary conditions (no communications) 32c 33 do j=1,nflow 34 p(:,j)=p(iper(:),j) 35 enddo 36c 37c slave has masters value, for abc we need to rotate it 38c (if this is a vector only no SCALARS) 39 if((iabc==1)) !are there any axisym bc's 40 & call rotabc(p(1,2), iBC, 'out') 41c 42c.... clear the vector 43c 44 q=zero 45c 46c.... Perform Matrix-Vector (AP) product 47c 48 do i = 1, nshg 49c 50 tmp1 = 0. 51 tmp2 = 0. 52 tmp3 = 0. 53 tmp4 = 0. 54 tmp5 = 0. 55c 56 do k = col(i), col(i+1)-1 57 j = row(k) 58c 59 tmp1 = tmp1 + lhsK(1 ,k)*p(j,1) 60 1 + lhsK(6 ,k)*p(j,2) 61 2 + lhsK(11,k)*p(j,3) 62 3 + lhsK(16,k)*p(j,4) 63 4 + lhsK(21,k)*p(j,5) 64 tmp2 = tmp2 + lhsK(2 ,k)*p(j,1) 65 1 + lhsK(7 ,k)*p(j,2) 66 2 + lhsK(12,k)*p(j,3) 67 3 + lhsK(17,k)*p(j,4) 68 4 + lhsK(22,k)*p(j,5) 69 tmp3 = tmp3 + lhsK(3 ,k)*p(j,1) 70 1 + lhsK(8 ,k)*p(j,2) 71 2 + lhsK(13,k)*p(j,3) 72 3 + lhsK(18,k)*p(j,4) 73 4 + lhsK(23,k)*p(j,5) 74 tmp4 = tmp4 + lhsK(4 ,k)*p(j,1) 75 1 + lhsK(9 ,k)*p(j,2) 76 2 + lhsK(14,k)*p(j,3) 77 3 + lhsK(19,k)*p(j,4) 78 4 + lhsK(24,k)*p(j,5) 79 tmp5 = tmp5 + lhsK(5 ,k)*p(j,1) 80 1 + lhsK(10,k)*p(j,2) 81 2 + lhsK(15,k)*p(j,3) 82 3 + lhsK(20,k)*p(j,4) 83 4 + lhsK(25,k)*p(j,5) 84c 85 enddo 86 87 q(i,1) = q(i,1) + tmp1 88 q(i,2) = q(i,2) + tmp2 89 q(i,3) = q(i,3) + tmp3 90 q(i,4) = q(i,4) + tmp4 91 q(i,5) = q(i,5) + tmp5 92 enddo 93 94 95 96c 97 p = q 98c 99c 100c.... --------------------> communications <------------------------- 101c 102c 103 if((iabc==1)) !are there any axisym bc's 104 & call rotabc(p(1,2), iBC, 'in ') 105c 106 if (numpe > 1) then 107c 108c.... send slave's copy of uBrg to the master 109c 110 call commu (p , ilwork, nflow , 'in ') 111c 112c.... nodes treated on another processor are eliminated 113c 114 numtask = ilwork(1) 115 itkbeg = 1 116 117 do itask = 1, numtask 118 119 iacc = ilwork (itkbeg + 2) 120 numseg = ilwork (itkbeg + 4) 121 122 if (iacc .eq. 0) then 123 do is = 1,numseg 124 isgbeg = ilwork (itkbeg + 3 + 2*is) 125 lenseg = ilwork (itkbeg + 4 + 2*is) 126 isgend = isgbeg + lenseg - 1 127 p(isgbeg:isgend,:) = zero 128 enddo 129 endif 130 131 itkbeg = itkbeg + 4 + 2*numseg 132 133 enddo 134 endif 135c 136 return 137 end 138 139