[430] | 1 | !The Full Polymorphic Package |
---|
| 2 | !Copyright (C) Etienne Forest and Frank Schmidt |
---|
| 3 | ! See file a_scratch_size |
---|
| 4 | |
---|
| 5 | module da_arrays |
---|
| 6 | use precision_constants |
---|
| 7 | use my_own_1D_TPSA |
---|
| 8 | use scratch_size |
---|
| 9 | implicit none |
---|
| 10 | public |
---|
| 11 | integer lda,lea,lia,lst |
---|
| 12 | integer,private,parameter::nmax=400 |
---|
| 13 | real(dp),private,parameter::tiny=1e-20_dp |
---|
| 14 | private ludcmp_nr,lubksb_nr |
---|
| 15 | ! johan |
---|
| 16 | ! integer,parameter::lno=2,lnv=6,lnomax=8,lnvmax=9,lstmax=800500,ldamax=16000,leamax=5000,liamax=50000 |
---|
| 17 | ! |
---|
| 18 | integer,parameter::lno=200,lnv=100,lnomax=8,lnvmax=9,lstmax=800500,ldamax=16000,leamax=5000,liamax=50000 |
---|
| 19 | logical(lp) :: reallocate = .true. |
---|
| 20 | logical(lp) :: notallocated = .true. |
---|
| 21 | logical(lp),parameter::etiennefix=.true. |
---|
| 22 | |
---|
| 23 | |
---|
| 24 | real(dp), allocatable,dimension(:)::cc ! lst |
---|
| 25 | integer, allocatable,dimension(:)::i_1,i_2 ! lst |
---|
| 26 | integer, allocatable,dimension(:)::ie1,ie2,ieo !lea |
---|
| 27 | integer, allocatable,dimension(:)::ia1,ia2 ! 0:lia |
---|
| 28 | integer, allocatable,dimension(:)::idano,idanv,idapo,idalm,idall ! lda |
---|
| 29 | character(10), allocatable,dimension(:)::daname ! daname(lda)*10 |
---|
| 30 | logical(lp), allocatable,dimension(:)::allvec !lda |
---|
| 31 | |
---|
| 32 | integer nda_dab |
---|
| 33 | integer :: ndamaxi=0 |
---|
| 34 | real(dp),TARGET :: total_da_size = 1.D38 !c_300 |
---|
| 35 | integer nst0,nomax,nvmax,nmmax,nocut,lfi |
---|
| 36 | real(dp) facint(0:lno) |
---|
| 37 | integer nhole |
---|
| 38 | integer,TARGET :: lda_used =15000 |
---|
| 39 | |
---|
| 40 | contains |
---|
| 41 | |
---|
| 42 | subroutine alloc_all(no,nv) |
---|
| 43 | implicit none |
---|
| 44 | integer no,nv |
---|
| 45 | if(reallocate) then |
---|
| 46 | call dealloc_all |
---|
| 47 | call danum0(no,nv) |
---|
| 48 | call alloc_ |
---|
| 49 | notallocated=.false. |
---|
| 50 | endif |
---|
| 51 | if(notallocated) then |
---|
| 52 | lia=liamax |
---|
| 53 | lst=lstmax |
---|
| 54 | lda=ldamax |
---|
| 55 | lea=leamax |
---|
| 56 | notallocated=.false. |
---|
| 57 | call alloc_ |
---|
| 58 | endif |
---|
| 59 | end subroutine alloc_all |
---|
| 60 | |
---|
| 61 | subroutine alloc_ |
---|
| 62 | implicit none |
---|
| 63 | integer i |
---|
| 64 | allocate(cc(lst)) |
---|
| 65 | allocate(i_1(lst));allocate(i_2(lst)); |
---|
| 66 | allocate(ie1(lea));allocate(ie2(lea));allocate(ieo(lea)); |
---|
| 67 | allocate(ia1(0:lia));allocate(ia2(0:lia)); |
---|
| 68 | allocate(idano(lda));allocate(idanv(lda));allocate(idapo(lda)); |
---|
| 69 | allocate(idalm(lda));allocate(idall(lda)); |
---|
| 70 | allocate(daname(lda)); |
---|
| 71 | allocate(allvec(lda)); |
---|
| 72 | do i=1,lst |
---|
| 73 | cc(I)=0.0_dp ! ADDED BY ETIENNE |
---|
| 74 | i_1(i)=0 |
---|
| 75 | i_2(i)=0 |
---|
| 76 | enddo |
---|
| 77 | do i=1,lea |
---|
| 78 | ie1(i)=0 |
---|
| 79 | ie2(i)=0 |
---|
| 80 | ieo(i)=0 |
---|
| 81 | enddo |
---|
| 82 | do i=0,lia |
---|
| 83 | ia1(i)=0 |
---|
| 84 | ia2(i)=0 |
---|
| 85 | enddo |
---|
| 86 | do i=1,lda |
---|
| 87 | idano(i)=0 |
---|
| 88 | idanv(i)=0 |
---|
| 89 | idapo(i)=0 |
---|
| 90 | idalm(i)=0 |
---|
| 91 | idall(i)=0 |
---|
| 92 | enddo |
---|
| 93 | end subroutine alloc_ |
---|
| 94 | |
---|
| 95 | |
---|
| 96 | subroutine dealloc_all |
---|
| 97 | implicit none |
---|
| 98 | integer error,ipause,mypauses |
---|
| 99 | IF (ALLOCATED(cc)) THEN |
---|
| 100 | DEALLOCATE (cc, STAT = error) |
---|
| 101 | IF(ERROR==0) THEN |
---|
| 102 | ELSE |
---|
| 103 | ipause=mypauses(100," cc ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 104 | ENDIF |
---|
| 105 | ENDIF |
---|
| 106 | IF (ALLOCATED(i_1)) THEN |
---|
| 107 | DEALLOCATE (i_1, STAT = error) |
---|
| 108 | IF(ERROR==0) THEN |
---|
| 109 | ELSE |
---|
| 110 | ipause=mypauses(101," i_1 ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 111 | ENDIF |
---|
| 112 | ENDIF |
---|
| 113 | IF (ALLOCATED(i_2)) THEN |
---|
| 114 | DEALLOCATE (i_2, STAT = error) |
---|
| 115 | IF(ERROR==0) THEN |
---|
| 116 | ELSE |
---|
| 117 | ipause=mypauses(102," i_2 ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 118 | ENDIF |
---|
| 119 | ENDIF |
---|
| 120 | IF (ALLOCATED(ie1)) THEN |
---|
| 121 | DEALLOCATE (ie1, STAT = error) |
---|
| 122 | IF(ERROR==0) THEN |
---|
| 123 | ELSE |
---|
| 124 | ipause=mypauses(103," ie1 ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 125 | ENDIF |
---|
| 126 | ENDIF |
---|
| 127 | IF (ALLOCATED(ie2)) THEN |
---|
| 128 | DEALLOCATE (ie2, STAT = error) |
---|
| 129 | IF(ERROR==0) THEN |
---|
| 130 | ELSE |
---|
| 131 | ipause=mypauses(104," ie2 ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 132 | ENDIF |
---|
| 133 | ENDIF |
---|
| 134 | IF (ALLOCATED(ieo)) THEN |
---|
| 135 | DEALLOCATE (ieo, STAT = error) |
---|
| 136 | IF(ERROR==0) THEN |
---|
| 137 | ELSE |
---|
| 138 | ipause=mypauses(105," ieo ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 139 | ENDIF |
---|
| 140 | ENDIF |
---|
| 141 | IF (ALLOCATED(ia1)) THEN |
---|
| 142 | DEALLOCATE (ia1, STAT = error) |
---|
| 143 | IF(ERROR==0) THEN |
---|
| 144 | ELSE |
---|
| 145 | ipause=mypauses(106," ia1 ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 146 | ENDIF |
---|
| 147 | ENDIF |
---|
| 148 | IF (ALLOCATED(ia2)) THEN |
---|
| 149 | DEALLOCATE (ia2, STAT = error) |
---|
| 150 | IF(ERROR==0) THEN |
---|
| 151 | ELSE |
---|
| 152 | ipause=mypauses(107," ia2 ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 153 | ENDIF |
---|
| 154 | ENDIF |
---|
| 155 | IF (ALLOCATED(idano)) THEN |
---|
| 156 | DEALLOCATE (idano, STAT = error) |
---|
| 157 | IF(ERROR==0) THEN |
---|
| 158 | ELSE |
---|
| 159 | ipause=mypauses(108," idano ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 160 | ENDIF |
---|
| 161 | ENDIF |
---|
| 162 | IF (ALLOCATED(idanv)) THEN |
---|
| 163 | DEALLOCATE (idanv, STAT = error) |
---|
| 164 | IF(ERROR==0) THEN |
---|
| 165 | ELSE |
---|
| 166 | ipause=mypauses(109," idanv ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 167 | ENDIF |
---|
| 168 | ENDIF |
---|
| 169 | IF (ALLOCATED(idapo)) THEN |
---|
| 170 | DEALLOCATE (idapo, STAT = error) |
---|
| 171 | IF(ERROR==0) THEN |
---|
| 172 | ELSE |
---|
| 173 | ipause=mypauses(110," idapo ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 174 | ENDIF |
---|
| 175 | ENDIF |
---|
| 176 | IF (ALLOCATED(idalm)) THEN |
---|
| 177 | DEALLOCATE (idalm, STAT = error) |
---|
| 178 | IF(ERROR==0) THEN |
---|
| 179 | ELSE |
---|
| 180 | ipause=mypauses(111," idalm ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 181 | ENDIF |
---|
| 182 | ENDIF |
---|
| 183 | IF (ALLOCATED(idall)) THEN |
---|
| 184 | DEALLOCATE (idall, STAT = error) |
---|
| 185 | IF(ERROR==0) THEN |
---|
| 186 | ELSE |
---|
| 187 | ipause=mypauses(112," idall ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 188 | ENDIF |
---|
| 189 | ENDIF |
---|
| 190 | IF (ALLOCATED(daname)) THEN |
---|
| 191 | DEALLOCATE (daname, STAT = error) |
---|
| 192 | IF(ERROR==0) THEN |
---|
| 193 | ELSE |
---|
| 194 | ipause=mypauses(112," daname ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 195 | ENDIF |
---|
| 196 | ENDIF |
---|
| 197 | IF (ALLOCATED(allvec)) THEN |
---|
| 198 | DEALLOCATE (allvec, STAT = error) |
---|
| 199 | IF(ERROR==0) THEN |
---|
| 200 | ELSE |
---|
| 201 | ipause=mypauses(112," allvec ARRAY not DEALLOCATED : PROBLEMS") |
---|
| 202 | ENDIF |
---|
| 203 | ENDIF |
---|
| 204 | end subroutine dealloc_all |
---|
| 205 | |
---|
| 206 | subroutine danum(no,nv,numda) |
---|
| 207 | implicit none |
---|
| 208 | integer i,mm,no,numda,nv |
---|
| 209 | ! ***************************** |
---|
| 210 | ! |
---|
| 211 | ! THIS SUBROUTINE COMPUTES THE NUMBER OF MONOMIALS OF |
---|
| 212 | ! ORDER NO AND NUMBER OF VARIABLES NV |
---|
| 213 | ! |
---|
| 214 | numda = 1 |
---|
| 215 | mm = max(nv,no) |
---|
| 216 | ! |
---|
| 217 | do i=1,min(nv,no) |
---|
| 218 | numda = (numda*(mm+i))/i |
---|
| 219 | enddo |
---|
| 220 | end subroutine danum |
---|
| 221 | |
---|
| 222 | subroutine danum0(no,nv) |
---|
| 223 | use scratch_size |
---|
| 224 | use file_handler |
---|
| 225 | implicit none |
---|
| 226 | integer i,mm,no,nv,ldamin ,nvt |
---|
| 227 | integer mf |
---|
| 228 | real(dp) size |
---|
| 229 | ! ***************************** |
---|
| 230 | ! |
---|
| 231 | ! THIS SUBROUTINE COMPUTES THE NUMBER OF MONOMIALS OF |
---|
| 232 | ! ORDER NO AND NUMBER OF VARIABLES NV |
---|
| 233 | ! |
---|
| 234 | lea = 1 |
---|
| 235 | if(etiennefix) then |
---|
| 236 | nvt=nv |
---|
| 237 | else |
---|
| 238 | if(nv==1) then ! still unknown feature of daini which must be taken into account |
---|
| 239 | nvt=2 ! perhaps related to Ying Wu modification in daini |
---|
| 240 | else |
---|
| 241 | nvt=nv |
---|
| 242 | endif |
---|
| 243 | endif |
---|
| 244 | mm = max(nvt,no) |
---|
| 245 | do i=1,min(nvt,no) |
---|
| 246 | lea = (lea*(mm+i))/i |
---|
| 247 | enddo |
---|
| 248 | ldamin=2+no + 2+ 4 !+ nd2t*(2+ndum) ! (2+ no)daini + varf +varc + assignmap(tpsalie.f90) |
---|
| 249 | |
---|
| 250 | |
---|
| 251 | ldamin=ldamin+2 ! 2 + sum of ndumuser in assign tpsa.f90 |
---|
| 252 | |
---|
| 253 | |
---|
| 254 | lia=(no+1)**((nv+mod(nv,2))/2) |
---|
| 255 | lda=lda_used+ldamin |
---|
| 256 | lst=lda*lea |
---|
| 257 | |
---|
| 258 | if(lingyun_yang) then |
---|
| 259 | size=REAL(lst,kind=DP)*(8.0_dp/1024.0_dp**2+2.0_dp*4.0_dp/1024.0_dp**2) |
---|
| 260 | else |
---|
| 261 | size=(3.0_dp*REAL(lea,kind=DP)+2.0_dp*(REAL(lia,kind=DP)+1.0_dp)+5.0_dp*REAL(lda,kind=DP))*4.0_dp/1024.0_dp**2 |
---|
| 262 | size=size+10.0_dp*REAL(lda,kind=DP)/1024.0_dp**2+REAL(lda,kind=DP)/8.0_dp/1024.0_dp**2 |
---|
| 263 | size=size+REAL(lst,kind=DP)*(8.0_dp/1024.0_dp**2+2.0_dp*4.0_dp/1024.0_dp**2) |
---|
| 264 | endif |
---|
| 265 | if(size>total_da_size.or.printdainfo) then |
---|
| 266 | w_p=0 |
---|
| 267 | w_p%nc=13 |
---|
| 268 | w_p%fc='(12(1X,A72,/),(1X,A72))' |
---|
| 269 | write(w_p%c(1),'(a10,1x,i4,1x,i4)') " no,nv = ",no,nv |
---|
| 270 | write(w_p%c(2),'(a10,1x,i8)') " LEA = ",lea |
---|
| 271 | write(w_p%c(3),'(a24,1x,i8)') " ldamin (with nd2=6) = ",ldamin |
---|
| 272 | write(w_p%c(4),'(a8,1x,i8)') " lia = ",lia |
---|
| 273 | write(w_p%c(5),'(a8,1x,i8)') " lda = ",lda |
---|
| 274 | write(w_p%c(6),'(a8,1x,i8)') " lst = ",lst |
---|
| 275 | write(w_p%c(7),'(a14,1x,i8)') " ndamaxi = ",ndamaxi |
---|
| 276 | write(w_p%c(8),'(a18,1x,g21.14)') " size in Mbytes = ",size |
---|
| 277 | write(w_p%c(9),'(a25,1x,g21.14)') " Total_da_size Allowed = ",Total_da_size |
---|
| 278 | w_p%c(10)=" " |
---|
| 279 | w_p%c(11)="************************************" |
---|
| 280 | w_p%c(12)="* Execution Continues Nevertheless *" |
---|
| 281 | w_p%c(13)="************************************" |
---|
| 282 | ! call !write_e(1000) |
---|
| 283 | call kanalnummer(mf) |
---|
| 284 | open(unit=mf,file='too_big_da.txt') |
---|
| 285 | write(mf,*) "no,nv = ",no,nv |
---|
| 286 | write(mf,*) "LEA = ",lea |
---|
| 287 | write(mf,*) "ldamin (with nd2=6) = ",ldamin |
---|
| 288 | write(mf,*) "lia = ",lia |
---|
| 289 | write(mf,*) "lda = ",lda |
---|
| 290 | write(mf,*) "lst = ",lst |
---|
| 291 | write(mf,*) "ndamaxi = ",ndamaxi |
---|
| 292 | write(mf,*) "size in Mbytes = ",size |
---|
| 293 | write(mf,*) "Total_da_size Allowed = ",Total_da_size |
---|
| 294 | write(6,*) "no,nv = ",no,nv |
---|
| 295 | write(6,*) "LEA = ",lea |
---|
| 296 | write(6,*) "ldamin (with nd2=6) = ",ldamin |
---|
| 297 | write(6,*) "lia = ",lia |
---|
| 298 | write(6,*) "lda = ",lda |
---|
| 299 | write(6,*) "lst = ",lst |
---|
| 300 | write(6,*) "ndamaxi = ",ndamaxi |
---|
| 301 | write(6,*) "size in Mbytes = ",size |
---|
| 302 | write(6,*) "Total_da_size Allowed = ",Total_da_size |
---|
| 303 | close(mf) |
---|
| 304 | endif |
---|
| 305 | |
---|
| 306 | |
---|
| 307 | end subroutine danum0 |
---|
| 308 | |
---|
| 309 | !!!!! some da stuff |
---|
| 310 | |
---|
| 311 | real(dp) function bran(xran) |
---|
| 312 | implicit none |
---|
| 313 | ! ************************************ |
---|
| 314 | ! |
---|
| 315 | ! VERY SIMPLE RANDOM NUMBER GENERATOR |
---|
| 316 | ! |
---|
| 317 | !----------------------------------------------------------------------------- |
---|
| 318 | ! |
---|
| 319 | real(dp) xran |
---|
| 320 | ! |
---|
| 321 | xran = xran + 10.0_dp |
---|
| 322 | if(xran.gt.1e4_dp) xran = xran - 9999.12345e0_dp |
---|
| 323 | bran = abs(sin(xran)) |
---|
| 324 | bran = 10*bran |
---|
| 325 | bran = bran - int(bran) |
---|
| 326 | ! IF(BRAN.LT. c_0_1) BRAN = BRAN + c_0_1 |
---|
| 327 | ! |
---|
| 328 | return |
---|
| 329 | end function bran |
---|
| 330 | |
---|
| 331 | subroutine matinv(a,ai,n,nmx,ier) |
---|
| 332 | |
---|
| 333 | implicit none |
---|
| 334 | |
---|
| 335 | integer i,ier,j,n,nmx |
---|
| 336 | integer,dimension(nmax)::indx |
---|
| 337 | real(dp) d |
---|
| 338 | real(dp),dimension(nmx,nmx)::a,ai |
---|
| 339 | real(dp),dimension(nmax,nmax)::aw |
---|
| 340 | ! |
---|
| 341 | ! if((.not.C_%STABLE_DA)) then |
---|
| 342 | ! if(c_%watch_user) then |
---|
| 343 | ! write(6,*) "big problem in dabnew ", sqrt(crash) |
---|
| 344 | ! endif |
---|
| 345 | ! return |
---|
| 346 | ! endif |
---|
| 347 | |
---|
| 348 | aw(1:n,1:n) = a(1:n,1:n) |
---|
| 349 | |
---|
| 350 | call ludcmp_nr(aw,n,nmax,indx,d,ier) |
---|
| 351 | if (ier .eq. 132) return |
---|
| 352 | |
---|
| 353 | ai(1:n,1:n) = 0.0_dp |
---|
| 354 | ! forall (i = 1:n) ai(i,i) = one |
---|
| 355 | do i=1,n |
---|
| 356 | ai(i,i) = 1.0_dp |
---|
| 357 | enddo |
---|
| 358 | |
---|
| 359 | do j=1,n |
---|
| 360 | call lubksb_nr(aw,n,nmax,indx,ai(1,j),nmx) |
---|
| 361 | enddo |
---|
| 362 | |
---|
| 363 | end subroutine matinv |
---|
| 364 | |
---|
| 365 | |
---|
| 366 | ! |
---|
| 367 | subroutine ludcmp_nr(a,n,np,indx,d,ier) |
---|
| 368 | implicit none |
---|
| 369 | ! ************************************ |
---|
| 370 | ! |
---|
| 371 | ! THIS SUBROUTINE DECOMPOSES A MATRIX INTO LU FORMAT |
---|
| 372 | ! INPUT A: NXN MATRIX - WILL BE OVERWRITTEN BY THE LU DECOMP. |
---|
| 373 | ! NP: PHYSICAL DIMENSION OF A |
---|
| 374 | ! INDX: ROW PERMUTATION VECTOR |
---|
| 375 | ! D: EVEN OR ODD ROW INTERCHANGES |
---|
| 376 | ! |
---|
| 377 | ! REFERENCE: NUMERICAL RECIPIES BY PRESS ET AL (CAMBRIDGE) PG. 35 |
---|
| 378 | ! |
---|
| 379 | !----------------------------------------------------------------------------- |
---|
| 380 | ! |
---|
| 381 | integer i,ier,imax,j,k,n,np |
---|
| 382 | integer,dimension(np)::indx |
---|
| 383 | real(dp) aamax,d,dum,sum |
---|
| 384 | real(dp),dimension(np,np)::a |
---|
| 385 | real(dp),dimension(nmax)::vv |
---|
| 386 | ! |
---|
| 387 | ! if((.not.C_%STABLE_DA)) then |
---|
| 388 | ! if(c_%watch_user) then |
---|
| 389 | ! write(6,*) "big problem in dabnew ", sqrt(crash) |
---|
| 390 | ! endif |
---|
| 391 | ! return |
---|
| 392 | ! endif |
---|
| 393 | ier=0 |
---|
| 394 | d=1.0_dp |
---|
| 395 | do i=1,n |
---|
| 396 | aamax=0.0_dp |
---|
| 397 | do j=1,n |
---|
| 398 | if(abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) |
---|
| 399 | enddo |
---|
| 400 | if(aamax.eq.0.0_dp) then |
---|
| 401 | ier=132 |
---|
| 402 | return |
---|
| 403 | endif |
---|
| 404 | vv(i)=1.0_dp/aamax |
---|
| 405 | enddo |
---|
| 406 | do j=1,n |
---|
| 407 | if(j.gt.1) then |
---|
| 408 | do i=1,j-1 |
---|
| 409 | sum=a(i,j) |
---|
| 410 | if(i.gt.1) then |
---|
| 411 | do k=1,i-1 |
---|
| 412 | sum=sum-a(i,k)*a(k,j) |
---|
| 413 | enddo |
---|
| 414 | a(i,j)=sum |
---|
| 415 | endif |
---|
| 416 | enddo |
---|
| 417 | endif |
---|
| 418 | aamax=0.0_dp |
---|
| 419 | do i=j,n |
---|
| 420 | sum=a(i,j) |
---|
| 421 | if (j.gt.1) then |
---|
| 422 | do k=1,j-1 |
---|
| 423 | sum=sum-a(i,k)*a(k,j) |
---|
| 424 | enddo |
---|
| 425 | a(i,j)=sum |
---|
| 426 | endif |
---|
| 427 | dum=vv(i)*abs(sum) |
---|
| 428 | if(dum.ge.aamax) then |
---|
| 429 | imax=i |
---|
| 430 | aamax=dum |
---|
| 431 | endif |
---|
| 432 | enddo |
---|
| 433 | if (j.ne.imax) then |
---|
| 434 | do k=1,n |
---|
| 435 | dum=a(imax,k) |
---|
| 436 | a(imax,k)=a(j,k) |
---|
| 437 | a(j,k)=dum |
---|
| 438 | enddo |
---|
| 439 | d=-d |
---|
| 440 | vv(imax)=vv(j) |
---|
| 441 | endif |
---|
| 442 | indx(j)=imax |
---|
| 443 | if(j.ne.n) then |
---|
| 444 | if(a(j,j).eq.0.0_dp) a(j,j)=tiny |
---|
| 445 | dum=1.0_dp/a(j,j) |
---|
| 446 | do i=j+1,n |
---|
| 447 | a(i,j)=a(i,j)*dum |
---|
| 448 | enddo |
---|
| 449 | endif |
---|
| 450 | enddo |
---|
| 451 | if(a(n,n).eq.0.0_dp) a(n,n)=tiny |
---|
| 452 | return |
---|
| 453 | end subroutine ludcmp_nr |
---|
| 454 | ! |
---|
| 455 | subroutine lubksb_nr(a,n,np,indx,b,nmx) |
---|
| 456 | implicit none |
---|
| 457 | ! ************************************ |
---|
| 458 | ! |
---|
| 459 | ! THIS SUBROUTINE SOLVES SET OF LINEAR EQUATIONS AX=B, |
---|
| 460 | ! INPUT A: NXN MATRIX IN lu FORM GIVEN BY ludcmp_nr |
---|
| 461 | ! NP: PHYSICAL DIMENSION OF A |
---|
| 462 | ! INDX: ROW PERMUTATION VECTOR |
---|
| 463 | ! D: EVEN OR ODD ROW INTERCHANGES |
---|
| 464 | ! B: RHS OF LINEAR EQUATION - WILL BE OVERWRITTEN BY X |
---|
| 465 | ! |
---|
| 466 | ! REFERENCE: NUMERICAL RECIPIES BY PRESS ET AL (CAMBRIDGE) PG. 36 |
---|
| 467 | ! |
---|
| 468 | !----------------------------------------------------------------------------- |
---|
| 469 | ! |
---|
| 470 | integer i,ii,j,ll,n,nmx,np |
---|
| 471 | integer,dimension(np)::indx |
---|
| 472 | real(dp) sum |
---|
| 473 | real(dp),dimension(np,np)::a |
---|
| 474 | real(dp),dimension(nmx)::b |
---|
| 475 | ! |
---|
| 476 | ! if((.not.C_%STABLE_DA)) then |
---|
| 477 | ! if(c_%watch_user) then |
---|
| 478 | ! write(6,*) "big problem in dabnew ", sqrt(crash) |
---|
| 479 | ! endif |
---|
| 480 | ! return |
---|
| 481 | ! endif |
---|
| 482 | ii = 0 |
---|
| 483 | do i=1,n |
---|
| 484 | ll = indx(i) |
---|
| 485 | sum = b(ll) |
---|
| 486 | b(ll) = b(i) |
---|
| 487 | if(ii.ne.0) then |
---|
| 488 | do j=ii,i-1 |
---|
| 489 | sum = sum-a(i,j)*b(j) |
---|
| 490 | enddo |
---|
| 491 | else if (sum.ne.0.0_dp) then |
---|
| 492 | ii = i |
---|
| 493 | endif |
---|
| 494 | b(i)=sum |
---|
| 495 | enddo |
---|
| 496 | do i=n,1,-1 |
---|
| 497 | sum=b(i) |
---|
| 498 | if(i.lt.n) then |
---|
| 499 | do j=i+1,n |
---|
| 500 | sum = sum-a(i,j)*b(j) |
---|
| 501 | enddo |
---|
| 502 | endif |
---|
| 503 | |
---|
| 504 | b(i)=sum/a(i,i) |
---|
| 505 | |
---|
| 506 | enddo |
---|
| 507 | return |
---|
| 508 | end subroutine lubksb_nr |
---|
| 509 | ! |
---|
| 510 | |
---|
| 511 | |
---|
| 512 | end module da_arrays |
---|