1 2!*********************************************************** 3! mtx_colm,mtx_rowp,mtx_mtx 4!********************************************************** 5 subroutine mtx_colm(a_colm,a_rowp, 6 & b_colm,b_rowp, 7 & c_colm,m,l,n) 8 integer,intent(in) :: m,l,n 9 integer,intent(in) :: a_colm(m+1),b_colm(l+1) 10 integer,intent(inout) :: c_colm(m+1) 11 integer,intent(in) :: a_rowp(a_colm(m+1)-1) 12 integer,intent(in) :: b_rowp(b_colm(l+1)-1) 13 14 integer :: i,j,k,p,q,mnnz 15 integer :: itemp(max(m,l)) 16 17 mnnz = 1 18 itemp = 0 19 do i=1,m 20 c_colm(i) = mnnz 21 do k=a_colm(i),a_colm(i+1)-1 22 j=a_rowp(k) 23 do p=b_colm(j),b_colm(j+1)-1 24 q=b_rowp(p) 25 if (itemp(q).ne.i) then 26 itemp(q)=i 27 mnnz = mnnz + 1 28 end if 29 enddo 30 enddo 31 enddo 32 c_colm(m+1) = mnnz 33 34 end subroutine!mtx_colm 35 36 subroutine mtx_rowp(a_colm,a_rowp,b_colm,b_rowp, 37 & c_colm,c_rowp,m,l,n,diagflag) 38 integer,intent(in) :: m,l,n 39 integer,intent(in) :: a_colm(m+1),b_colm(l+1) 40 integer,intent(in) :: c_colm(m+1) 41 integer,intent(inout) :: c_rowp(c_colm(m+1)-1) 42 integer,intent(in) :: a_rowp(a_colm(m+1)-1) 43 integer,intent(in) :: b_rowp(b_colm(l+1)-1) 44 logical,intent(in) :: diagflag 45 46 integer :: i,j,k,p,q,mnnz 47 integer :: itemp(max(m,l)) 48 49 mnnz = 0 50 itemp = 0 51 do i=1,m 52 if (diagflag) then 53 mnnz = mnnz + 1 54 c_rowp(mnnz) = i 55 itemp(i)=i 56 endif 57 do k=a_colm(i),a_colm(i+1)-1 58 j=a_rowp(k) 59 do p=b_colm(j),b_colm(j+1)-1 60 q=b_rowp(p) 61 if (itemp(q).ne.i) then 62 itemp(q)=i 63 mnnz = mnnz + 1 64 c_rowp(mnnz) = q 65 end if 66 enddo 67 enddo 68 enddo 69 70 end subroutine !mtx_rowp 71 72 subroutine mtx_mtx(a_colm,a_rowp,a_mtx, 73 & b_colm,b_rowp,b_mtx, 74 & c_colm,c_rowp,c_mtx, 75 & m,l,n) 76 integer,intent(in) :: m,l,n 77 integer,intent(in) :: a_colm(m+1),b_colm(l+1) 78 integer,intent(in) :: c_colm(m+1) 79 integer,intent(in) :: c_rowp(c_colm(m+1)-1) 80 integer,intent(in) :: a_rowp(a_colm(m+1)-1) 81 integer,intent(in) :: b_rowp(b_colm(l+1)-1) 82 83 real(kind=8),intent(in) :: a_mtx(a_colm(m+1)-1) 84 real(kind=8),intent(in) :: b_mtx(b_colm(l+1)-1) 85 real(kind=8),intent(inout) :: c_mtx(c_colm(m+1)-1) 86 87 integer :: i,j,k,p,q 88 integer :: itemp(max(m,l)) 89 real(kind=8) :: rtemp 90 91 c_mtx = 0 92 do i=1,m 93 do k=c_colm(i),c_colm(i+1)-1 94 itemp(c_rowp(k)) = k 95 enddo 96 do k=a_colm(i),a_colm(i+1)-1 97 j = a_rowp(k) 98 rtemp = a_mtx(k) 99 do p = b_colm(j),b_colm(j+1)-1 100 q = b_rowp(p) 101 c_mtx(itemp(q)) = c_mtx(itemp(q)) + rtemp*b_mtx(p) 102 enddo 103 enddo 104 enddo 105 106 end subroutine !mtx_mtx 107 108!************************************************************ 109! ramg_calcIAI 110! do matrix multiplication for I^T A I 111! actual operation is: A2_ij = I^T_ip A_pq I_qj 112! where p,q are repeated indicies 113! 114!************************************************************ 115 subroutine ramg_calcITAI(level1,level2,maxstopsign) 116 use ramg_data 117 implicit none 118 119 integer,intent(in) :: level1,level2 120 logical,intent(inout) :: maxstopsign 121 122 integer :: i,j,k,m,p,q,r,n,s,p0 123 integer,dimension(amg_nshg(level1)) :: itemp 124 integer,dimension(amg_nshg(level2)) :: itemp2 125 integer :: mnnz 126 127 integer,dimension(amg_nshg(level2)+1) :: ITA_colm 128 integer,dimension(:),allocatable :: ITA_rowp 129 real(kind=8),dimension(:),allocatable :: ITA 130 131 integer :: mem_err, mem_err_s 132 133 real(kind=8) :: rtp 134 character :: fname*80 135 real,dimension(10) :: cpusec 136 137 call cpu_time(cpusec(1)) 138 139 IF (level1.eq.0) THEN 140 ! reduced serial test here, DO NOT USE IT unless RECOVER 141 ! from extract. 142 call mtx_colm(I_fc_colm(level1)%p,I_fc_rowp(level1)%p, 143! & amg_A_colm(level1)%p,amg_A_rowp(level1)%p, 144 & reducecolm,reducerowp, 145 & ITA_colm,amg_nshg(level2),amg_nshg(level1), 146 & amg_nshg(level1)) 147 mnnz = ITA_colm(amg_nshg(level2)+1)-1 148 allocate(ITA_rowp(mnnz),stat=mem_err) 149 allocate(ITA(mnnz),stat=mem_err) 150 ! IT * A 151 ! rowp 152 call mtx_rowp(I_fc_colm(level1)%p,I_fc_rowp(level1)%p, 153! & amg_A_colm(level1)%p,amg_A_rowp(level1)%p, 154 & reducecolm,reducerowp, 155 & ITA_colm,ITA_rowp, 156 & amg_nshg(level2),amg_nshg(level1), 157 & amg_nshg(level1),.false.) 158 ! IT * A 159 ! value 160 call mtx_mtx(I_fc_colm(level1)%p,I_fc_rowp(level1)%p, 161 & I_fc(level1)%p, 162! & amg_A_colm(level1)%p,amg_A_rowp(level1)%p, 163! & amg_A_lhs(level1)%p, 164 & reducecolm,reducerowp,reducelhs, 165 & ITA_colm,ITA_rowp,ITA, 166 & amg_nshg(level2),amg_nshg(level1), 167 & amg_nshg(level1)) 168 call cpu_time(cpusec(2)) 169 170 ELSE 171 172 ! IT * A 173 ! colm and nnz 174 call mtx_colm(I_fc_colm(level1)%p,I_fc_rowp(level1)%p, 175 & amg_A_colm(level1)%p,amg_A_rowp(level1)%p, 176 & ITA_colm,amg_nshg(level2),amg_nshg(level1), 177 & amg_nshg(level1)) 178 mnnz = ITA_colm(amg_nshg(level2)+1)-1 179 allocate(ITA_rowp(mnnz),stat=mem_err) 180 allocate(ITA(mnnz),stat=mem_err) 181 ! IT * A 182 ! rowp 183 call mtx_rowp(I_fc_colm(level1)%p,I_fc_rowp(level1)%p, 184 & amg_A_colm(level1)%p,amg_A_rowp(level1)%p, 185 & ITA_colm,ITA_rowp, 186 & amg_nshg(level2),amg_nshg(level1), 187 & amg_nshg(level1),.false.) 188 ! IT * A 189 ! value 190 call mtx_mtx(I_fc_colm(level1)%p,I_fc_rowp(level1)%p, 191 & I_fc(level1)%p, 192 & amg_A_colm(level1)%p,amg_A_rowp(level1)%p, 193 & amg_A_lhs(level1)%p, 194 & ITA_colm,ITA_rowp,ITA, 195 & amg_nshg(level2),amg_nshg(level1), 196 & amg_nshg(level1)) 197 call cpu_time(cpusec(2)) 198 199 ENDIF 200 201 ! ( IT*A ) * I 202 ! colm and nnz 203 call ramg_allocate(level2,amg_nshg(level2),0,1) 204 call mtx_colm(ITA_colm,ITA_rowp,I_cf_colm(level1)%p, 205 & I_cf_rowp(level1)%p,amg_A_colm(level2)%p, 206 & amg_nshg(level2),amg_nshg(level1), 207 & amg_nshg(level2)) 208 209 call cpu_time(cpusec(3)) 210 mnnz = amg_A_colm(level2)%p(amg_nshg(level2)+1)-1 211 call ramg_allocate(level2,0,mnnz,1) 212 ! ( IT*A ) * I 213 ! rowp 214 call mtx_rowp(ITA_colm,ITA_rowp,I_cf_colm(level1)%p, 215 & I_cf_rowp(level1)%p,amg_A_colm(level2)%p, 216 & amg_A_rowp(level2)%p, 217 & amg_nshg(level2),amg_nshg(level1), 218 & amg_nshg(level2),.true.) 219 call cpu_time(cpusec(4)) 220 ! ( IT*A ) * I 221 ! value 222 call mtx_mtx(ITA_colm,ITA_rowp,ITA,I_cf_colm(level1)%p, 223 & I_cf_rowp(level1)%p,I_cf(level1)%p, 224 & amg_A_colm(level2)%p, 225 & amg_A_rowp(level2)%p,amg_A_lhs(level2)%p, 226 & amg_nshg(level2),amg_nshg(level1), 227 & amg_nshg(level2)) 228 229 deallocate(ITA_rowp,stat=mem_err) 230 deallocate(ITA,stat=mem_err) 231 232 call cpu_time(cpusec(5)) 233 234 amg_nshg(level2+1) = 0 235 236 ! put on right hand side 237 call ramg_calcIvFC(level1,level2,amg_A_rhs(level1)%p, 238 & amg_A_rhs(level2)%p) 239 240 call cpu_time(cpusec(10)) 241 242 do i=1,amg_nshg(level2) 243 amg_ppeDiag(level2)%p(i) = 244 & 1.0/sqrt(amg_A_lhs(level2)%p(amg_A_colm(level2)%p(i),1)) 245 enddo 246 do i=1,amg_nshg(level2) 247 do m=amg_A_colm(level2)%p(i),amg_A_colm(level2)%p(i+1)-1 248 j = amg_A_rowp(level2)%p(m) 249 amg_A_lhs(level2)%p(m,1)=amg_A_lhs(level2)%p(m,1)* 250 & amg_ppeDiag(level2)%p(i)*amg_ppeDiag(level2)%p(j) 251 enddo 252 enddo 253 254! do i=1,amg_nshg(level1) 255! do j=I_cf_colm(level1)%p(i),I_cf_colm(level1)%p(i+1)-1 256! k = I_cf_rowp(level1)%p(j) 257! I_cf(level1)%p(j) = 258! & I_cf(level1)%p(j)*amg_ppeDiag(level2)%p(k) 259! enddo 260! enddo 261 262! do i=1,amg_nshg(level2) 263! do j=I_fc_colm(level1)%p(i),I_fc_colm(level1)%p(i+1)-1 264! I_cf(level1)%p(j)= 265! & I_cf(level1)%p(j)*amg_ppeDiag(level2)%p(i) 266! enddo 267! enddo 268 269 end subroutine ! ramg_calcIAI 270 271 272