1!$Id:$ 2 subroutine cassem(ad, au, al, s, ir, jc, ld, nst, alfl, 3 & bycol,diagin,all) 4 5! * * F E A P * * A Finite Element Analysis Program 6 7!.... Copyright (c) 1984-2017: Regents of the University of California 8! All rights reserved 9 10!-----[--.----+----.----+----.-----------------------------------------] 11! Purpose: Assemble a sparse matrix A into a compressed array 12! 'ad', 'au', 'al' 13 14! Inputs: 15! s(*,*) - Element matrix to assemble 16! ir(*) - Location of non-zero entries in A by profile col/rows 17! jc(*) - Pointer array to find entries for equations 18! ld(*) - Local/global array to map 's' into 'A'. 19! nst - Size of 's' and 'ld' 20! alfl - Flag: if true assemble both upper and lower parts 21! if false assemble symmetric array in upper part 22! bycol - Sparse storage scheme by column if .true. 23! diagin - Includes diagonal in sparse assembly if .true. 24! all - Stores all terms in row/column if .true. 25 26! Outputs: 27! ad(*) - Diagonal part of A 28! au(*) - Sparse store of upper part of A 29! al(*) - Sparse store of lower part of A 30!-----[--.----+----.----+----.-----------------------------------------] 31 implicit none 32 33 include 'cdata.h' 34 include 'constant.h' 35 36 logical alfl,bycol,diagin,all 37 integer nst, i,j,k,n 38 integer inz,inza, jc(*),ir(*), ld(*) 39 real*8 ad(*), au(*), al(*), s(nst,nst) 40 41 save 42 43! Compact assembly of profile matrix 44 45 do i = 1,nst 46 n = ld(i) 47 if( n.ge.1) then 48 49! Assemble total array 50 51 if(all) then 52 53! Assemble by columns 54 55 if(bycol) then 56 do j = 1,nst 57 k = ld(j) 58 if(k.gt.0 .and. n.gt.1) then 59 inz = inza(jc(n-1)+1, jc(n), ir, k, n) 60 ad(inz) = ad(inz) + s(j,i) 61 elseif(k.gt.0 .and. n.eq.1) then 62 inz = inza( one, jc(n), ir, k, n) 63 ad(inz) = ad(inz) + s(j,i) 64 endif 65 end do ! j 66 67! Assemble by rows (no reordering) 68 69 else 70 do j = 1,nst 71 k = ld(j) 72 if(k.gt.0 .and. n.gt.1) then 73 inz = inza(jc(n-1)+1, jc(n), ir, k, n) 74 ad(inz) = ad(inz) + s(i,j) 75 elseif(k.gt.0 .and. n.eq.1) then 76 inz = inza( one, jc(n), ir, k, n) 77 ad(inz) = ad(inz) + s(i,j) 78 endif 79 end do ! j 80 81 endif 82 83! Assemble upper/lower parts by columns 84 85 elseif(bycol) then 86 do j = 1,nst 87 k = ld(j) 88! Assemble including diagonal 89 if(diagin) then 90 if(k.gt.0 .and. k.le.n ) then 91 if(n.eq.1) then 92 ad(1) = ad(1) + s(j,i) 93 if(alfl) al(1) = al(1) + s(i,j) 94 else 95 inz = inza(jc(n-1)+1, jc(n), ir, k, n) 96 ad(inz) = ad(inz) + s(j,i) 97 if(alfl) al(inz) = al(inz) + s(i,j) 98 endif 99 endif 100! Assemble excluding diagonal 101 else 102 if(k.gt.0 .and. k.lt.n ) then 103 inz = inza(jc(n-1)+1, jc(n), ir, k, n) 104 au(inz) = au(inz) + s(j,i) 105 if(alfl) al(inz) = al(inz) + s(i,j) 106 endif 107 endif 108 end do ! j 109 110! Assemble upper/lower parts by rows 111 112 else 113 do j = 1,nst 114 k = ld(j) 115 if(k.gt.0) then 116! Assemble including diagonal 117 if(diagin) then 118 if( k.ge.n ) then 119 inz = inza(jc(n+neq), jc(n+neq+1)-1, ir, k, n) 120 ad(inz) = ad(inz) + s(j,i) 121 if(alfl) al(inz) = al(inz) + s(i,j) 122 endif 123! Assemble excluding diagonal 124 else 125 if( k.gt.n ) then 126 inz = inza(jc(n+neq), jc(n+neq+1)-1, ir, k, n) 127 au(inz) = au(inz) + s(j,i) 128 if(alfl) al(inz) = al(inz) + s(i,j) 129 endif 130 endif 131 endif 132 end do ! j 133 endif 134 135! Assemble diagonal for .not.diagin cases 136 137 if(.not.diagin) then 138 do j = 1,nst 139 if(ld(j).eq.n) ad(n) = ad(n) + s(i,j) 140 end do ! j 141 endif 142 143 endif 144 145 end do ! i 146 147 end 148