xref: /phasta/phSolver/common/fillsparse.f (revision 712d3df0b59ebebaaeaea358162c8d2c043c6e08)
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