1!======================================================================! 2! ! 3!======================================================================! 4SUBROUTINE D6DOT(T,A,B,N) 5 IMPLICIT NONE 6 7 INTEGER :: jj 8 INTEGER :: l 9 INTEGER :: N 10 DOUBLE PRECISION :: T(9) 11 DOUBLE PRECISION :: A(9,*) 12 DOUBLE PRECISION :: B(9,*) 13 !---------------------------------------------------------------------- 14 ! 15 ! spdot1 performs inner product of sparse vectors 16 ! 17 ! 18 ! #coded by t.arakawa of RIST on 040510 19 ! 20 !---------------------------------------------------------------------- 21 DO l = 1, 9 22 T(l) = 0.0D0 23 ENDDO 24 DO jj = 1, N 25 T(1) = T(1) + A(1,jj)*B(1,jj) + A(4,jj)*B(4,jj) + A(7,jj)*B(7,jj) 26 T(2) = T(2) + A(2,jj)*B(1,jj) + A(5,jj)*B(4,jj) + A(8,jj)*B(7,jj) 27 T(3) = T(3) + A(3,jj)*B(1,jj) + A(6,jj)*B(4,jj) + A(9,jj)*B(7,jj) 28 T(4) = T(4) + A(1,jj)*B(2,jj) + A(4,jj)*B(5,jj) + A(7,jj)*B(8,jj) 29 T(5) = T(5) + A(2,jj)*B(2,jj) + A(5,jj)*B(5,jj) + A(8,jj)*B(8,jj) 30 T(6) = T(6) + A(3,jj)*B(2,jj) + A(6,jj)*B(5,jj) + A(9,jj)*B(8,jj) 31 T(7) = T(7) + A(1,jj)*B(3,jj) + A(4,jj)*B(6,jj) + A(7,jj)*B(9,jj) 32 T(8) = T(8) + A(2,jj)*B(3,jj) + A(5,jj)*B(6,jj) + A(8,jj)*B(9,jj) 33 T(9) = T(9) + A(3,jj)*B(3,jj) + A(6,jj)*B(6,jj) + A(9,jj)*B(9,jj) 34 ENDDO 35END SUBROUTINE D6DOT 36!======================================================================! 37! ! 38!======================================================================! 39SUBROUTINE D6DOTL(T,A,B,N) 40 IMPLICIT NONE 41 42 INTEGER :: jj 43 INTEGER :: l 44 INTEGER :: N 45 DOUBLE PRECISION :: T(6) 46 DOUBLE PRECISION :: A(9,*) 47 DOUBLE PRECISION :: B(9,*) 48 !---------------------------------------------------------------------- 49 ! 50 ! spdot1 performs inner product of sparse vectors 51 ! 52 ! 53 ! #coded by t.arakawa of RIST on 040510 54 ! 55 !---------------------------------------------------------------------- 56 !$dir max_trips(6) 57 DO l = 1, 6 58 T(l) = 0.0D0 59 ENDDO 60 DO jj = 1, N 61 T(1) = T(1) + A(1,jj)*B(1,jj) + A(4,jj)*B(4,jj) + A(7,jj)*B(7,jj) 62 T(2) = T(2) + A(2,jj)*B(1,jj) + A(5,jj)*B(4,jj) + A(8,jj)*B(7,jj) 63 T(3) = T(3) + A(2,jj)*B(2,jj) + A(5,jj)*B(5,jj) + A(8,jj)*B(8,jj) 64 T(4) = T(4) + A(3,jj)*B(1,jj) + A(6,jj)*B(4,jj) + A(9,jj)*B(7,jj) 65 T(5) = T(5) + A(3,jj)*B(2,jj) + A(6,jj)*B(5,jj) + A(9,jj)*B(8,jj) 66 T(6) = T(6) + A(3,jj)*B(3,jj) + A(6,jj)*B(6,jj) + A(9,jj)*B(9,jj) 67 ENDDO 68END SUBROUTINE D6DOTL 69!======================================================================! 70! ! 71!======================================================================! 72SUBROUTINE D6SDOT(Wi,A,B,N) 73 IMPLICIT NONE 74 75 INTEGER :: jj 76 INTEGER :: N 77 DOUBLE PRECISION :: Wi(3) 78 DOUBLE PRECISION :: A(3,*) 79 DOUBLE PRECISION :: B(9,*) 80 !---------------------------------------------------------------------- 81 ! 82 ! spdot1 performs inner product of sparse vectors 83 ! 84 ! 85 ! #coded by t.arakawa of RIST on 040510 86 ! 87 !---------------------------------------------------------------------- 88 DO jj = 1, N 89 Wi(1) = Wi(1) - A(1,jj)*B(1,jj) - A(2,jj)*B(4,jj) - A(3,jj)*B(7,jj) 90 Wi(2) = Wi(2) - A(1,jj)*B(2,jj) - A(2,jj)*B(5,jj) - A(3,jj)*B(8,jj) 91 Wi(3) = Wi(3) - A(1,jj)*B(3,jj) - A(2,jj)*B(6,jj) - A(3,jj)*B(9,jj) 92 ENDDO 93END SUBROUTINE D6SDOT 94!======================================================================! 95! ! 96!======================================================================! 97SUBROUTINE IDNTTY(Neqns,Invp,Iperm) 98 IMPLICIT NONE 99 100 INTEGER :: i 101 INTEGER :: IDBg1 102 INTEGER :: Neqns 103 INTEGER :: Invp(*) 104 INTEGER :: Iperm(*) 105 COMMON /DEBUG / IDBg1 106 107 i = 1 108 DO WHILE ( i<=Neqns ) 109 WRITE (6,*) 'invp(', i, ')' 110 READ (5,*) Invp(i) 111 IF ( Invp(i)==0 ) THEN 112 DO i = 1, Neqns 113 Invp(i) = i 114 Iperm(i) = i 115 ENDDO 116 RETURN 117 ELSEIF ( Invp(i)<0 ) THEN 118 READ (11,*) (Invp(i),i=1,Neqns) 119 DO i = 1, Neqns 120 Iperm(Invp(i)) = i 121 ENDDO 122 GOTO 99999 123 ELSE 124 i = i + 1 125 ENDIF 126 ENDDO 127 DO i = 1, Neqns 128 Iperm(Invp(i)) = i 129 ENDDO 130 RETURN 13199999 END SUBROUTINE IDNTTY 132 !======================================================================! 133 ! ! 134 !======================================================================! 135SUBROUTINE NUSOL6(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop) 136 IMPLICIT NONE 137 138 INTEGER :: i 139 INTEGER :: j 140 INTEGER :: joc 141 INTEGER :: k 142 INTEGER :: ke 143 INTEGER :: ks 144 INTEGER :: Neqns 145 INTEGER :: Nstop 146 INTEGER :: Xlnzr(*) 147 INTEGER :: Colno(*) 148 INTEGER :: Iperm(*) 149 !GP: DEBUG 13May04 wk(3 ---> wk(6, b(3 ---> b(6, diag(6 ---> diag(21, 150 !GP: DEBUG 13May04 zln(9 ---> zln(36, dsln(9 ---> dsln(36 151 ! double precision zln(9,*),diag(6,*),b(3,*),wk(3,*),dsln(9,*) 152 DOUBLE PRECISION :: Zln(36,*) 153 DOUBLE PRECISION :: Diag(21,*) 154 DOUBLE PRECISION :: B(6,*) 155 DOUBLE PRECISION :: Wk(6,*) 156 DOUBLE PRECISION :: Dsln(36,*) 157 ! forward 158 DO i = 1, Neqns 159 Wk(1,i) = B(1,Iperm(i)) 160 Wk(2,i) = B(2,Iperm(i)) 161 Wk(3,i) = B(3,Iperm(i)) 162 Wk(4,i) = B(4,Iperm(i)) 163 Wk(5,i) = B(5,Iperm(i)) 164 Wk(6,i) = B(6,Iperm(i)) 165 ENDDO 166 joc = 1 167 DO i = 1, Neqns 168 ks = Xlnzr(i) 169 ke = Xlnzr(i+1) - 1 170 IF ( ke>=ks ) CALL S6PDOT(Wk(1,i),Wk,Zln,Colno,ks,ke) 171 IF ( i>Nstop ) THEN 172 ! call d6sdot(wk(1,i),wk(1,nstop),dsln(1,joc),i-nstop) 173 CALL DXSDOT(6,Wk(1,i),Wk(1,Nstop),Dsln(1,joc),i-Nstop) 174 joc = joc + i - Nstop 175 ENDIF 176 ENDDO 177 DO i = 1, Neqns 178 Wk(2,i) = Wk(2,i) - Wk(1,i)*Diag(2,i) 179 Wk(3,i) = Wk(3,i) - Wk(1,i)*Diag(4,i) - Wk(2,i)*Diag(5,i) 180 Wk(4,i) = Wk(4,i) - Wk(1,i)*Diag(7,i) - Wk(2,i)*Diag(8,i) - Wk(3,i)*Diag(9,i) 181 Wk(5,i) = Wk(5,i) - Wk(1,i)*Diag(11,i) - Wk(2,i)*Diag(12,i) - Wk(3,i)*Diag(13,i) - Wk(4,i)*Diag(14,i) 182 Wk(6,i) = Wk(6,i) - Wk(1,i)*Diag(16,i) - Wk(2,i)*Diag(17,i)& 183 - Wk(3,i)*Diag(18,i) - Wk(4,i)*Diag(19,i) - Wk(6,i)& 184 *Diag(20,i) 185 Wk(1,i) = Wk(1,i)*Diag(1,i) 186 Wk(2,i) = Wk(2,i)*Diag(3,i) 187 Wk(3,i) = Wk(3,i)*Diag(6,i) 188 Wk(4,i) = Wk(4,i)*Diag(10,i) 189 Wk(5,i) = Wk(5,i)*Diag(15,i) 190 Wk(6,i) = Wk(6,i)*Diag(21,i) 191 Wk(5,i) = Wk(5,i) - Wk(6,i)*Diag(20,i) 192 Wk(4,i) = Wk(4,i) - Wk(6,i)*Diag(19,i) - Wk(5,i)*Diag(14,i) 193 Wk(3,i) = Wk(3,i) - Wk(6,i)*Diag(18,i) - Wk(5,i)*Diag(13,i) - Wk(4,i)*Diag(9,i) 194 Wk(2,i) = Wk(2,i) - Wk(6,i)*Diag(17,i) - Wk(5,i)*Diag(12,i) - Wk(4,i)*Diag(8,i) - Wk(3,i)*Diag(5,i) 195 Wk(1,i) = Wk(1,i) - Wk(6,i)*Diag(16,i) - Wk(5,i)*Diag(11,i)& 196 - Wk(4,i)*Diag(7,i) - Wk(3,i)*Diag(4,i) - Wk(2,i)& 197 *Diag(2,i) 198 ENDDO 199 ! back ward 200 DO i = Neqns, 1, -1 201 IF ( i>=Nstop ) THEN 202 DO j = i - 1, Nstop, -1 203 joc = joc - 1 204 Wk(1,j) = Wk(1,j) - Wk(1,i)*Dsln(1,joc) - Wk(2,i)& 205 *Dsln(2,joc) - Wk(3,i)*Dsln(3,joc) - Wk(4,i)& 206 *Dsln(4,joc) - Wk(5,i)*Dsln(5,joc) - Wk(6,i)& 207 *Dsln(6,joc) 208 Wk(2,j) = Wk(2,j) - Wk(1,i)*Dsln(7,joc) - Wk(2,i)& 209 *Dsln(8,joc) - Wk(3,i)*Dsln(9,joc) - Wk(4,i)& 210 *Dsln(10,joc) - Wk(5,i)*Dsln(11,joc) - Wk(6,i)& 211 *Dsln(12,joc) 212 Wk(3,j) = Wk(3,j) - Wk(1,i)*Dsln(13,joc) - Wk(2,i)& 213 *Dsln(14,joc) - Wk(3,i)*Dsln(15,joc) - Wk(4,i)& 214 *Dsln(16,joc) - Wk(5,i)*Dsln(17,joc) - Wk(6,i)& 215 *Dsln(18,joc) 216 Wk(4,j) = Wk(4,j) - Wk(1,i)*Dsln(19,joc) - Wk(2,i)& 217 *Dsln(20,joc) - Wk(3,i)*Dsln(21,joc) - Wk(4,i)& 218 *Dsln(22,joc) - Wk(5,i)*Dsln(23,joc) - Wk(6,i)& 219 *Dsln(24,joc) 220 Wk(5,j) = Wk(5,j) - Wk(1,i)*Dsln(25,joc) - Wk(2,i)& 221 *Dsln(26,joc) - Wk(3,i)*Dsln(27,joc) - Wk(4,i)& 222 *Dsln(28,joc) - Wk(5,i)*Dsln(29,joc) - Wk(6,i)& 223 *Dsln(30,joc) 224 Wk(6,j) = Wk(6,j) - Wk(1,i)*Dsln(31,joc) - Wk(2,i)& 225 *Dsln(32,joc) - Wk(3,i)*Dsln(33,joc) - Wk(4,i)& 226 *Dsln(34,joc) - Wk(5,i)*Dsln(35,joc) - Wk(6,i)& 227 *Dsln(36,joc) 228 ENDDO 229 ENDIF 230 ks = Xlnzr(i) 231 ke = Xlnzr(i+1) - 1 232 IF ( ke>=ks ) THEN 233 DO k = ks, ke 234 j = Colno(k) 235 Wk(1,j) = Wk(1,j) - Wk(1,i)*Zln(1,joc) - Wk(2,i)& 236 *Zln(2,joc) - Wk(3,i)*Zln(3,joc) - Wk(4,i)& 237 *Zln(4,joc) - Wk(5,i)*Zln(5,joc) - Wk(6,i)& 238 *Zln(6,joc) 239 Wk(2,j) = Wk(2,j) - Wk(1,i)*Zln(7,joc) - Wk(2,i)& 240 *Zln(8,joc) - Wk(3,i)*Zln(9,joc) - Wk(4,i)& 241 *Zln(10,joc) - Wk(5,i)*Zln(11,joc) - Wk(6,i)& 242 *Zln(12,joc) 243 Wk(3,j) = Wk(3,j) - Wk(1,i)*Zln(13,joc) - Wk(2,i)& 244 *Zln(14,joc) - Wk(3,i)*Zln(15,joc) - Wk(4,i)& 245 *Zln(16,joc) - Wk(5,i)*Zln(17,joc) - Wk(6,i)& 246 *Zln(18,joc) 247 Wk(4,j) = Wk(4,j) - Wk(1,i)*Zln(19,joc) - Wk(2,i)& 248 *Zln(20,joc) - Wk(3,i)*Zln(21,joc) - Wk(4,i)& 249 *Zln(22,joc) - Wk(5,i)*Zln(23,joc) - Wk(6,i)& 250 *Zln(24,joc) 251 Wk(5,j) = Wk(5,j) - Wk(1,i)*Zln(25,joc) - Wk(2,i)& 252 *Zln(26,joc) - Wk(3,i)*Zln(27,joc) - Wk(4,i)& 253 *Zln(28,joc) - Wk(5,i)*Zln(29,joc) - Wk(6,i)& 254 *Zln(30,joc) 255 Wk(6,j) = Wk(6,j) - Wk(1,i)*Zln(31,joc) - Wk(2,i)& 256 *Zln(32,joc) - Wk(3,i)*Zln(33,joc) - Wk(4,i)& 257 *Zln(34,joc) - Wk(5,i)*Zln(35,joc) - Wk(6,i)& 258 *Zln(36,joc) 259 ENDDO 260 ENDIF 261 ENDDO 262 ! permutaion 263 DO i = 1, Neqns 264 B(1,Iperm(i)) = Wk(1,i) 265 B(2,Iperm(i)) = Wk(2,i) 266 B(3,Iperm(i)) = Wk(3,i) 267 B(4,Iperm(i)) = Wk(4,i) 268 B(5,Iperm(i)) = Wk(5,i) 269 B(6,Iperm(i)) = Wk(6,i) 270 ENDDO 271END SUBROUTINE NUSOL6 272!======================================================================! 273! ! 274!======================================================================! 275SUBROUTINE PRT(Ip,N) 276 IMPLICIT NONE 277 278 INTEGER :: i 279 INTEGER :: Ip 280 INTEGER :: N 281 DIMENSION Ip(N) 282 WRITE (6,99001) (Ip(i),i=1,N) 28399001 FORMAT (10(2x,i4)) 284END SUBROUTINE PRT 285 286SUBROUTINE VLCPY1(A,N) 287 IMPLICIT NONE 288 289 DOUBLE PRECISION :: A 290 INTEGER :: N 291 DIMENSION A(N) 292 INTEGER :: i 293 INTEGER :: j 294 A(N) = 0 295END SUBROUTINE VLCPY1 296 297!======================================================================! 298! ! 299!======================================================================! 300SUBROUTINE VERIF0(Neqns,Ndeg,Nttbr,Irow,Jcol,Val,Rhs,X) 301 IMPLICIT NONE 302 303 DOUBLE PRECISION :: err 304 DOUBLE PRECISION :: rel 305 DOUBLE PRECISION :: Rhs 306 DOUBLE PRECISION :: Val 307 DOUBLE PRECISION :: X 308 INTEGER :: i 309 INTEGER :: Irow 310 INTEGER :: j 311 INTEGER :: Jcol 312 INTEGER :: k 313 INTEGER :: l 314 INTEGER :: m 315 INTEGER :: Ndeg 316 INTEGER :: Neqns 317 INTEGER :: Nttbr 318 DIMENSION Irow(*), Jcol(*), Val(Ndeg,Ndeg,*), Rhs(Ndeg,*), X(Ndeg,*) 319 !---------------------------------------------------------------------- 320 ! 321 ! verify the solution(symmetric matrix) 322 ! 323 !---------------------------------------------------------------------- 324 rel = 0.0D0 325 DO i = 1, Neqns 326 DO l = 1, Ndeg 327 rel = rel + DABS(Rhs(l,i)) 328 ENDDO 329 ENDDO 330 DO k = 1, Nttbr 331 i = Irow(k) 332 j = Jcol(k) 333 DO l = 1, Ndeg 334 DO m = 1, Ndeg 335 Rhs(l,i) = Rhs(l,i) - Val(l,m,k)*X(m,j) 336 IF ( i/=j ) Rhs(l,j) = Rhs(l,j) - Val(m,l,k)*X(m,i) 337 ENDDO 338 ENDDO 339 ENDDO 340 err = 0.0D0 341 DO i = 1, Neqns 342 DO l = 1, Ndeg 343 err = err + DABS(Rhs(l,i)) 344 ENDDO 345 ENDDO 346 WRITE (6,99001) err, rel, err/rel 347 !WINDEBUG 348 ! write(16,6000) err,rel,err/rel 34999001 FORMAT (' ***verification***(symmetric)'/& 350 'norm(Ax-b) = ',& 351 1pd20.10/'norm(b) = ',& 352 1pd20.10/'norm(Ax-b)/norm(b) = ',1pd20.10) 353END SUBROUTINE VERIF0 354 355!======================================================================! 356! ! 357!======================================================================! 358SUBROUTINE V6PROD(Zln,Diag,Zz,N) 359 IMPLICIT NONE 360 361 DOUBLE PRECISION :: Diag 362 DOUBLE PRECISION :: Zln 363 DOUBLE PRECISION :: Zz 364 INTEGER :: i 365 INTEGER :: N 366 DIMENSION Zln(9,N), Diag(6,N), Zz(9,N) 367 DO i = 1, N 368 Zz(4,i) = Zln(4,i) - Zln(1,i)*Diag(2,i) 369 Zz(7,i) = Zln(7,i) - Zln(1,i)*Diag(4,i) - Zz(4,i)*Diag(5,i) 370 Zz(1,i) = Zln(1,i)*Diag(1,i) 371 Zz(4,i) = Zz(4,i)*Diag(3,i) 372 Zz(7,i) = Zz(7,i)*Diag(6,i) 373 Zz(4,i) = Zz(4,i) - Zz(7,i)*Diag(5,i) 374 Zz(1,i) = Zz(1,i) - Zz(4,i)*Diag(2,i) - Zz(7,i)*Diag(4,i) 375 ! 376 Zz(5,i) = Zln(5,i) - Zln(2,i)*Diag(2,i) 377 Zz(8,i) = Zln(8,i) - Zln(2,i)*Diag(4,i) - Zz(5,i)*Diag(5,i) 378 Zz(2,i) = Zln(2,i)*Diag(1,i) 379 Zz(5,i) = Zz(5,i)*Diag(3,i) 380 Zz(8,i) = Zz(8,i)*Diag(6,i) 381 Zz(5,i) = Zz(5,i) - Zz(8,i)*Diag(5,i) 382 Zz(2,i) = Zz(2,i) - Zz(5,i)*Diag(2,i) - Zz(8,i)*Diag(4,i) 383 ! 384 Zz(6,i) = Zln(6,i) - Zln(3,i)*Diag(2,i) 385 Zz(9,i) = Zln(9,i) - Zln(3,i)*Diag(4,i) - Zz(6,i)*Diag(5,i) 386 Zz(3,i) = Zln(3,i)*Diag(1,i) 387 Zz(6,i) = Zz(6,i)*Diag(3,i) 388 Zz(9,i) = Zz(9,i)*Diag(6,i) 389 Zz(6,i) = Zz(6,i) - Zz(9,i)*Diag(5,i) 390 Zz(3,i) = Zz(3,i) - Zz(6,i)*Diag(2,i) - Zz(9,i)*Diag(4,i) 391 ENDDO 392END SUBROUTINE V6PROD 393