159599516SKenneth E. Jansen subroutine fillsparseI( iens, xKebe, lhsK, 259599516SKenneth E. Jansen & xGoC, lhsP, 359599516SKenneth E. Jansen 1 row, col) 459599516SKenneth E. Jansenc 559599516SKenneth E. Jansenc 659599516SKenneth E. Jansenc 759599516SKenneth E. Jansen include "common.h" 859599516SKenneth E. Jansen real*8 xKebe(npro,9,nshl,nshl), xGoC(npro,4,nshl,nshl) 959599516SKenneth E. Jansen integer ien(npro,nshl), col(nshg+1), row(nshg*nnz) 1059599516SKenneth E. Jansen real*8 lhsK(9,nnz_tot), lhsP(4,nnz_tot) 1159599516SKenneth E. Jansenc 1259599516SKenneth E. Jansen integer aa, b, c, e, i, k, n 1359599516SKenneth E. Jansenc 1459599516SKenneth E. Jansen integer sparseloc 1559599516SKenneth E. Jansen 1659599516SKenneth E. Jansen integer iens(npro,nshl) 1759599516SKenneth E. Jansenc 1859599516SKenneth E. Jansenc prefer to show explicit absolute value needed for cubic modes and 1959599516SKenneth E. Jansenc higher rather than inline abs on pointer as in past versions 2059599516SKenneth E. Jansenc iens is the signed ien array ien is unsigned 2159599516SKenneth E. Jansenc 2259599516SKenneth E. Jansen ien=abs(iens) 2359599516SKenneth E. Jansenc 2459599516SKenneth E. Jansenc.... Accumulate the lhs 2559599516SKenneth E. Jansenc 2659599516SKenneth E. Jansen do e = 1, npro ! loop over the elements 2759599516SKenneth E. Jansen do aa = 1, nshl ! loop over the local equation numbers 2859599516SKenneth E. Jansen i = ien(e,aa) ! finds the global equation number or 2959599516SKenneth E. Jansen ! block-row of our matrix 3059599516SKenneth E. Jansen c = col(i) ! starting point to look for the matching column 3159599516SKenneth E. Jansen n = col(i+1) - c !length of the list of entries in rowp 3259599516SKenneth E. Jansen do b = 1, nshl ! local variable number tangent respect 3359599516SKenneth E. Jansen ! to 3459599516SKenneth E. Jansenc function that searches row until it finds the match that gives the 3559599516SKenneth E. Jansenc global equation number 3659599516SKenneth E. Jansen 3759599516SKenneth E. Jansen k = sparseloc( row(c), n, ien(e,b) ) + c-1 3859599516SKenneth E. Jansenc 3959599516SKenneth E. Jansenc * * 4059599516SKenneth E. Jansenc dimension egmass(npro,ndof,nenl,ndof,nenl) 4159599516SKenneth E. Jansenc 4259599516SKenneth E. Jansenc compressible lhsT(1:5,1:5,k)=lhsT(1:5,1:5,k)+egmass(e,1:5,aa,1:5,b) 4359599516SKenneth E. Jansenc 4459599516SKenneth E. Jansen lhsK(1,k) = lhsK(1,k) + xKebe(e,1,aa,b) 4559599516SKenneth E. Jansen lhsK(2,k) = lhsK(2,k) + xKebe(e,2,aa,b) 4659599516SKenneth E. Jansen lhsK(3,k) = lhsK(3,k) + xKebe(e,3,aa,b) 4759599516SKenneth E. Jansen lhsK(4,k) = lhsK(4,k) + xKebe(e,4,aa,b) 4859599516SKenneth E. Jansen lhsK(5,k) = lhsK(5,k) + xKebe(e,5,aa,b) 4959599516SKenneth E. Jansen lhsK(6,k) = lhsK(6,k) + xKebe(e,6,aa,b) 5059599516SKenneth E. Jansen lhsK(7,k) = lhsK(7,k) + xKebe(e,7,aa,b) 5159599516SKenneth E. Jansen lhsK(8,k) = lhsK(8,k) + xKebe(e,8,aa,b) 5259599516SKenneth E. Jansen lhsK(9,k) = lhsK(9,k) + xKebe(e,9,aa,b) 5359599516SKenneth E. Jansenc 5459599516SKenneth E. Jansen lhsP(1,k) = lhsP(1,k) + xGoC(e,1,aa,b) 5559599516SKenneth E. Jansen lhsP(2,k) = lhsP(2,k) + xGoC(e,2,aa,b) 5659599516SKenneth E. Jansen lhsP(3,k) = lhsP(3,k) + xGoC(e,3,aa,b) 5759599516SKenneth E. Jansen lhsP(4,k) = lhsP(4,k) + xGoC(e,4,aa,b) 5859599516SKenneth E. Jansen enddo 5959599516SKenneth E. Jansen enddo 6059599516SKenneth E. Jansen enddo 6159599516SKenneth E. Jansenc 6259599516SKenneth E. Jansenc.... end 6359599516SKenneth E. Jansenc 6459599516SKenneth E. Jansen return 6559599516SKenneth E. Jansen end 6659599516SKenneth E. Jansen 6759599516SKenneth E. Jansen 6859599516SKenneth E. Jansen subroutine fillsparseC( iens, EGmass, lhsK, 6959599516SKenneth E. Jansen 1 row, col) 7059599516SKenneth E. Jansenc 7159599516SKenneth E. Jansenc----------------------------------------------------------- 7259599516SKenneth E. Jansenc This routine fills up the spasely stored LHS mass matrix 7359599516SKenneth E. Jansenc 7459599516SKenneth E. Jansenc Nahid Razmra, Spring 2000. (Sparse Matrix) 7559599516SKenneth E. Jansenc----------------------------------------------------------- 7659599516SKenneth E. Jansenc 7759599516SKenneth E. Jansenc 7859599516SKenneth E. Jansen 7959599516SKenneth E. Jansen include "common.h" 8059599516SKenneth E. Jansen 8159599516SKenneth E. Jansen real*8 EGmass(npro,nedof,nedof) 8259599516SKenneth E. Jansen integer ien(npro,nshl), col(nshg+1), row(nnz*nshg) 8359599516SKenneth E. Jansen real*8 lhsK(nflow*nflow,nnz_tot) 8459599516SKenneth E. Jansen 8559599516SKenneth E. Jansenc 8659599516SKenneth E. Jansen integer aa, b, c, e, i, k, n, n1 8759599516SKenneth E. Jansen integer f, g, r, s, t 8859599516SKenneth E. Jansenc 8959599516SKenneth E. Jansen integer sparseloc 9059599516SKenneth E. Jansen 9159599516SKenneth E. Jansen integer iens(npro,nshl) 9259599516SKenneth E. Jansenc 9359599516SKenneth E. Jansenc prefer to show explicit absolute value needed for cubic modes and 9459599516SKenneth E. Jansenc higher rather than inline abs on pointer as in past versions 9559599516SKenneth E. Jansenc iens is the signed ien array ien is unsigned 9659599516SKenneth E. Jansenc 9759599516SKenneth E. Jansen ien=abs(iens) 9859599516SKenneth E. Jansenc 9959599516SKenneth E. Jansenc.... Accumulate the lhsK 10059599516SKenneth E. Jansenc 10159599516SKenneth E. Jansen do e = 1, npro 102*513954efSKenneth E. Jansen do aa = 1, nshl !loop over matrix block column index 103*513954efSKenneth E. Jansen i = ien(e,aa) !get mapping from aath node of element e to 104*513954efSKenneth E. Jansen ! the ith row of lhsk 105*513954efSKenneth E. Jansen 106*513954efSKenneth E. Jansen c = col(i) !Get the mapping from the ith row of lhsk to 107*513954efSKenneth E. Jansen ! to the corresponding location in row 108*513954efSKenneth E. Jansen n = col(i+1) - c !number of nonzero blocks in the row 109*513954efSKenneth E. Jansen r = (aa-1)*nflow !starting index of the ath node in EGmass 11059599516SKenneth E. Jansen do b = 1, nshl 111*513954efSKenneth E. Jansen s = (b-1)*nflow !starting index of the bth node's 112*513954efSKenneth E. Jansen !contribution to node aa. 11359599516SKenneth E. Jansen k = sparseloc( row(c), n, ien(e,b) ) + c-1 114*513954efSKenneth E. Jansen !Find the index of row which corresponds to node b of 115*513954efSKenneth E. Jansen !element e. This is where contributions from EGmass 116*513954efSKenneth E. Jansen !will actually get stored. 117*513954efSKenneth E. Jansen 118*513954efSKenneth E. Jansen do g = 1, nflow !loop over columns in a block 11959599516SKenneth E. Jansen t = (g-1)*nflow 120*513954efSKenneth E. Jansen do f = 1, nflow !loop over rows in a block 12159599516SKenneth E. Jansen lhsK(t+f,k) = lhsK(t+f,k) + EGmass(e,r+f,s+g) 122*513954efSKenneth E. Jansen enddo 123*513954efSKenneth E. Jansen enddo 124*513954efSKenneth E. Jansen enddo 125*513954efSKenneth E. Jansen enddo 126*513954efSKenneth E. Jansen enddo 12759599516SKenneth E. Jansen 12859599516SKenneth E. Jansen return 12959599516SKenneth E. Jansen end 13059599516SKenneth E. Jansen 131*513954efSKenneth E. Jansen subroutine fillsparseC_BC( iens, EGmass, 132*513954efSKenneth E. Jansen & lhsk, row, col) 133*513954efSKenneth E. Jansen! 134*513954efSKenneth E. Jansen!----------------------------------------------------------- 135*513954efSKenneth E. Jansen! This routine adds contributions from heat flux BCs to the 136*513954efSKenneth E. Jansen! spasely stored LHS mass matrix. This routine is modified 137*513954efSKenneth E. Jansen! from fillsparseC. 138*513954efSKenneth E. Jansen! 139*513954efSKenneth E. Jansen! Nicholas Mati, Summer 2014. (Sparse Matrix) 140*513954efSKenneth E. Jansen!----------------------------------------------------------- 141*513954efSKenneth E. Jansen! 142*513954efSKenneth E. Jansen include "common.h" 143*513954efSKenneth E. Jansen 144*513954efSKenneth E. Jansen real*8 EGmass(npro, nshl, nshl) !only contains term (5,5) 145*513954efSKenneth E. Jansen integer ien(npro,nshl), col(nshg+1), row(nnz*nshg) 146*513954efSKenneth E. Jansen real*8 lhsK(nflow*nflow,nnz_tot) 147*513954efSKenneth E. Jansen 148*513954efSKenneth E. Jansen integer aa, b, c, e, i, k, n, n1 149*513954efSKenneth E. Jansen integer f, g, r, s, t 150*513954efSKenneth E. Jansen 151*513954efSKenneth E. Jansen integer sparseloc 152*513954efSKenneth E. Jansen integer iens(npro,nshl) 153*513954efSKenneth E. Jansenc 154*513954efSKenneth E. Jansenc prefer to show explicit absolute value needed for cubic modes and 155*513954efSKenneth E. Jansenc higher rather than inline abs on pointer as in past versions 156*513954efSKenneth E. Jansenc iens is the signed ien array ien is unsigned 157*513954efSKenneth E. Jansenc 158*513954efSKenneth E. Jansen ien=abs(iens) 159*513954efSKenneth E. Jansen 160*513954efSKenneth E. Jansen !Accumulate the lhsK 161*513954efSKenneth E. Jansen do e = 1, npro 162*513954efSKenneth E. Jansen do aa = 1, nshl !loop over matrix block column index 163*513954efSKenneth E. Jansen i = ien(e,aa) !get mapping from aath node of element e to 164*513954efSKenneth E. Jansen ! the ith row of lhsk 165*513954efSKenneth E. Jansen 166*513954efSKenneth E. Jansen c = col(i) !Get the mapping from the ith row of lhsk to 167*513954efSKenneth E. Jansen ! to the corresponding location in row 168*513954efSKenneth E. Jansen n = col(i+1) - c !number of nonzero blocks in the row 169*513954efSKenneth E. Jansen! r = (aa-1)*nflow !starting index of the ath node in EGmass 170*513954efSKenneth E. Jansen do b = 1, nshl 171*513954efSKenneth E. Jansen! s = (b-1)*nflow !starting index of the bth node's 172*513954efSKenneth E. Jansen !contribution to node aa. 173*513954efSKenneth E. Jansen k = sparseloc( row(c), n, ien(e,b) ) + c-1 174*513954efSKenneth E. Jansen !Find the index of row which corresponds to node b of 175*513954efSKenneth E. Jansen !element e. This is where contributions from EGmass 176*513954efSKenneth E. Jansen !will actually get stored. 177*513954efSKenneth E. Jansen 178*513954efSKenneth E. Jansen lhsk(25, k) = lhsk(25, k) + EGmass(e, aa, b) 179*513954efSKenneth E. Jansen! do g = 1, nflow !loop over columns in a block 180*513954efSKenneth E. Jansen! t = (g-1)*nflow 181*513954efSKenneth E. Jansen! do f = 1, nflow !loop over rows in a block 182*513954efSKenneth E. Jansen! lhsK(t+f,k) = lhsK(t+f,k) + EGmass(e,r+f,s+g) 183*513954efSKenneth E. Jansen! enddo 184*513954efSKenneth E. Jansen! enddo 185*513954efSKenneth E. Jansen enddo !loop over node b 186*513954efSKenneth E. Jansen enddo !loop over node a 187*513954efSKenneth E. Jansen enddo !loop over elements 188*513954efSKenneth E. Jansen 189*513954efSKenneth E. Jansen !end 190*513954efSKenneth E. Jansen return 191*513954efSKenneth E. Jansen end 192*513954efSKenneth E. Jansen 193*513954efSKenneth E. Jansen 19459599516SKenneth E. Jansen subroutine fillsparseSclr( iens, xSebe, lhsS, 19559599516SKenneth E. Jansen 1 row, col) 196*513954efSKenneth E. Jansen 19759599516SKenneth E. Jansen include "common.h" 19859599516SKenneth E. Jansen real*8 xSebe(npro,nshl,nshl) 19959599516SKenneth E. Jansen integer ien(npro,nshl), col(nshg+1), row(nshg*nnz) 20059599516SKenneth E. Jansen real*8 lhsS(nnz_tot) 201*513954efSKenneth E. Jansen 20259599516SKenneth E. Jansen integer aa, b, c, e, i, k, n 203*513954efSKenneth E. Jansen 20459599516SKenneth E. Jansen integer sparseloc 20559599516SKenneth E. Jansen 20659599516SKenneth E. Jansen integer iens(npro,nshl) 20759599516SKenneth E. Jansenc 20859599516SKenneth E. Jansenc prefer to show explicit absolute value needed for cubic modes and 20959599516SKenneth E. Jansenc higher rather than inline abs on pointer as in past versions 21059599516SKenneth E. Jansenc iens is the signed ien array ien is unsigned 21159599516SKenneth E. Jansenc 21259599516SKenneth E. Jansen ien=abs(iens) 21359599516SKenneth E. Jansenc 21459599516SKenneth E. Jansenc.... Accumulate the lhs 21559599516SKenneth E. Jansenc 21659599516SKenneth E. Jansen do e = 1, npro 21759599516SKenneth E. Jansen do aa = 1, nshl 21859599516SKenneth E. Jansen i = ien(e,aa) 21959599516SKenneth E. Jansen c = col(i) 22059599516SKenneth E. Jansen n = col(i+1) - c 22159599516SKenneth E. Jansen do b = 1, nshl 22259599516SKenneth E. Jansen k = sparseloc( row(c), n, ien(e,b) ) + c-1 22359599516SKenneth E. Jansenc 22459599516SKenneth E. Jansen lhsS(k) = lhsS(k) + xSebe(e,aa,b) 22559599516SKenneth E. Jansen enddo 22659599516SKenneth E. Jansen enddo 22759599516SKenneth E. Jansen enddo 22859599516SKenneth E. Jansenc 22959599516SKenneth E. Jansenc.... end 23059599516SKenneth E. Jansenc 23159599516SKenneth E. Jansen return 23259599516SKenneth E. Jansen end 23359599516SKenneth E. Jansen 23459599516SKenneth E. Jansen integer function sparseloc( list, n, target ) 23559599516SKenneth E. Jansen 23659599516SKenneth E. Jansenc----------------------------------------------------------- 23759599516SKenneth E. Jansenc This function finds the location of the non-zero elements 23859599516SKenneth E. Jansenc of the LHS matrix in the sparsely stored matrix 23959599516SKenneth E. Jansenc lhsK(nflow*nflow,nnz*numnp) 24059599516SKenneth E. Jansenc 24159599516SKenneth E. Jansenc Nahid Razmara, Spring 2000. (Sparse Matrix) 24259599516SKenneth E. Jansenc----------------------------------------------------------- 24359599516SKenneth E. Jansen 24459599516SKenneth E. Jansen integer list(n), n, target 24559599516SKenneth E. Jansen integer rowvl, rowvh, rowv 24659599516SKenneth E. Jansen 24759599516SKenneth E. Jansenc 24859599516SKenneth E. Jansenc.... Initialize 24959599516SKenneth E. Jansenc 25059599516SKenneth E. Jansen rowvl = 1 25159599516SKenneth E. Jansen rowvh = n + 1 25259599516SKenneth E. Jansenc 25359599516SKenneth E. Jansenc.... do a binary search 25459599516SKenneth E. Jansenc 25559599516SKenneth E. Jansen100 if ( rowvh-rowvl .gt. 1 ) then 25659599516SKenneth E. Jansen rowv = ( rowvh + rowvl ) / 2 25759599516SKenneth E. Jansen if ( list(rowv) .gt. target ) then 25859599516SKenneth E. Jansen rowvh = rowv 25959599516SKenneth E. Jansen else 26059599516SKenneth E. Jansen rowvl = rowv 26159599516SKenneth E. Jansen endif 26259599516SKenneth E. Jansen goto 100 26359599516SKenneth E. Jansen endif 26459599516SKenneth E. Jansenc 26559599516SKenneth E. Jansenc.... return 26659599516SKenneth E. Jansenc 26759599516SKenneth E. Jansen sparseloc = rowvl 26859599516SKenneth E. Jansenc 26959599516SKenneth E. Jansen return 27059599516SKenneth E. Jansen end 27159599516SKenneth E. Jansen 272