Millepede-II V04-17-04
Dbandmatrix.f90
Go to the documentation of this file.
1
2! Code converted using TO_F90 by Alan Miller
3! Date: 2012-04-18 Time: 19:56:08
4
232
233
234! (1) Symmetric matrix W: decomposition, solution, inverse
235
244
245SUBROUTINE dchdec(w,n, aux)
246 USE mpdef
247
248 implicit none
249 INTEGER(mpi) :: i
250 INTEGER(mpi) :: ii
251 INTEGER(mpi) :: j
252 INTEGER(mpi) :: jj
253 INTEGER(mpi) :: k
254 INTEGER(mpi) :: kk
255 REAL(mpd) :: ratio
256
257 INTEGER(mpi), INTENT(IN) :: n
258 REAL(mpd), INTENT(IN OUT) :: w((n*n+n)/2)
259 REAL(mpd), INTENT(OUT) :: aux(n)
260
261 ! ...
262 DO i=1,n
263 aux(i)=16.0_mpd*w((i*i+i)/2) ! save diagonal elements
264 END DO
265 ii=0
266 DO i=1,n
267 ii=ii+i
268 IF(w(ii)+aux(i) /= aux(i)) THEN ! GT
269 w(ii)=1.0_mpd/w(ii) ! (I,I) div !
270 ELSE
271 w(ii)=0.0_mpd
272 END IF
273 jj=ii
274 DO j=i+1,n
275 ratio=w(i+jj)*w(ii) ! (I,J) (I,I)
276 kk=jj
277 DO k=j,n
278 w(kk+j)=w(kk+j)-w(kk+i)*ratio ! (K,J) (K,I)
279 kk=kk+k
280 END DO ! K
281 w(i+jj)=ratio ! (I,J)
282 jj=jj+j
283 END DO ! J
284 END DO ! I
285 RETURN
286END SUBROUTINE dchdec
287
294
295SUBROUTINE dchslv(w,n,b, x)
296 USE mpdef
297
298 implicit none
299 INTEGER(mpi) :: i
300 INTEGER(mpi) :: ii
301 INTEGER(mpi) :: k
302 INTEGER(mpi) :: kk
303 REAL(mpd) :: suma
304
305 INTEGER(mpi), INTENT(IN) :: n
306 REAL(mpd), INTENT(IN) :: w((n*n+n)/2)
307 REAL(mpd), INTENT(IN) :: b(n)
308 REAL(mpd), INTENT(OUT) :: x(n)
309
310 WRITE(*,*) 'before copy'
311 DO i=1,n
312 x(i)=b(i)
313 END DO
314 WRITE(*,*) 'after copy'
315 ii=0
316 DO i=1,n
317 suma=x(i)
318 DO k=1,i-1
319 suma=suma-w(k+ii)*x(k) ! (K,I)
320 END DO
321 x(i)=suma
322 ii=ii+i
323 END DO
324 WRITE(*,*) 'after forward'
325 DO i=n,1,-1
326 suma=x(i)*w(ii) ! (I,I)
327 kk=ii
328 DO k=i+1,n
329 suma=suma-w(kk+i)*x(k) ! (K,I)
330 kk=kk+k
331 END DO
332 x(i)=suma
333 ii=ii-i
334 END DO
335 WRITE(*,*) 'after backward'
336 RETURN
337END SUBROUTINE dchslv
338
344
345SUBROUTINE dchinv(w,n,v)
346 USE mpdef
347
348 implicit none
349 INTEGER(mpi) :: i
350 INTEGER(mpi) :: ii
351 INTEGER(mpi) :: j
352 INTEGER(mpi) :: k
353 INTEGER(mpi) :: l
354 INTEGER(mpi) :: m
355 REAL(mpd) :: suma
356
357 INTEGER(mpi), INTENT(IN) :: n
358 REAL(mpd), INTENT(IN) :: w((n*n+n)/2)
359 REAL(mpd), INTENT(OUT) :: v((n*n+n)/2)
360
361 ii=(n*n-n)/2
362 DO i=n,1,-1
363 suma=w(ii+i) ! (I,I)
364 DO j=i,1,-1
365 DO k=j+1,n
366 l=min(i,k)
367 m=max(i,k)
368 suma=suma-w(j+(k*k-k)/2)*v(l+(m*m-m)/2) ! (J,K) (I,K)
369 END DO
370 v(ii+j)=suma ! (I,J)
371 suma=0.0_mpd
372 END DO
373 ii=ii-i+1
374 END DO
375END SUBROUTINE dchinv
376
383
384REAL(mps) FUNCTION condes(w,n,aux)
385 USE mpdef
386
387 implicit none
388 REAL(mps) :: cond
389 INTEGER(mpi) :: i
390 INTEGER(mpi) :: ir
391 INTEGER(mpi) :: is
392 INTEGER(mpi) :: k
393 INTEGER(mpi) :: kk
394 REAL(mps) :: xla1
395 REAL(mps) :: xlan
396
397
398 INTEGER(mpi), INTENT(IN) :: n
399 REAL(mpd), INTENT(IN) :: w((n*n+n)/2)
400 REAL(mpd), INTENT(IN OUT) :: aux(n)
401
402 REAL(mpd) :: suma,sumu,sums
403
404 ir=1
405 is=1
406 DO i=1,n
407 IF(w((i*i+i)/2) < w((is*is+is)/2)) is=i ! largest Dii
408 IF(w((i*i+i)/2) > w((ir*ir+ir)/2)) ir=i ! smallest Dii
409 END DO
410
411 sumu=0.0 ! find smallest eigenvalue
412 sums=0.0
413 DO i=n,1,-1 ! backward
414 suma=0.0
415 IF(i == ir) suma=1.0_mpd
416 kk=(i*i+i)/2
417 DO k=i+1,n
418 suma=suma-w(kk+i)*aux(k) ! (K,I)
419 kk=kk+k
420 END DO
421 aux(i)=suma
422 sumu=sumu+aux(i)*aux(i)
423 END DO
424 xlan=real(w((ir*ir+ir)/2)*sqrt(sumu),mps)
425 IF(xlan /= 0.0) xlan=1.0/xlan
426
427 DO i=1,n
428 IF(i == is) THEN
429 sums=1.0_mpd
430 ELSE IF(i > is) THEN
431 sums=sums+w(is+(i*i-i)/2)**2
432 END IF
433 END DO ! is Ws
434 xla1=0.0
435 IF(w((is*is+is)/2) /= 0.0) xla1=real(sqrt(sums)/w((is*is+is)/2),mps)
436
437 cond=0.0
438 IF(xla1 > 0.0.AND.xlan > 0.0) cond=xla1/xlan
439 ! estimated condition
440 condes=cond
441END FUNCTION condes
442
443
444! (2) Symmetric band matrix: decomposition, solution, complete
445! inverse and band part of the inverse
459SUBROUTINE dbcdec(w,mp1,n, aux)
460 USE mpdef
461
462 implicit none
463 INTEGER(mpi) :: i
464 INTEGER(mpi) :: j
465 INTEGER(mpi) :: k
466 REAL(mpd) :: rxw
467 ! M=MP1-1 N*M(M-1) dot operations
468
469 INTEGER(mpi), INTENT(IN) :: n
470 INTEGER(mpi), INTENT(IN) :: mp1
471 REAL(mpd), INTENT(IN OUT) :: w(mp1,n)
472 REAL(mpd), INTENT(OUT) :: aux(n)
473 ! decompos
474 ! ...
475 DO i=1,n
476 aux(i)=16.0_mpd*w(1,i) ! save diagonal elements
477 END DO
478 DO i=1,n
479 IF(w(1,i)+aux(i) /= aux(i)) THEN
480 w(1,i)=1.0/w(1,i)
481 ELSE
482 w(1,i)=0.0_mpd ! singular
483 END IF
484 DO j=1,min(mp1-1,n-i)
485 rxw=w(j+1,i)*w(1,i)
486 DO k=1,min(mp1-1,n-i)+1-j
487 w(k,i+j)=w(k,i+j)-w(k+j,i)*rxw
488 END DO
489 w(j+1,i)=rxw
490 END DO
491 END DO
492 RETURN
493END SUBROUTINE dbcdec
494
502
503SUBROUTINE dbcslv(w,mp1,n,b, x)
504 USE mpdef
505
506 implicit none
507 INTEGER(mpi) :: i
508 INTEGER(mpi) :: j
509 REAL(mpd) :: rxw
510
511 INTEGER(mpi), INTENT(IN) :: n
512 INTEGER(mpi), INTENT(IN) :: mp1
513 REAL(mpd), INTENT(IN) :: w(mp1,n)
514 REAL(mpd), INTENT(IN) :: b(n)
515 REAL(mpd), INTENT(OUT) :: x(n)
516 ! N*(2M-1) dot operations
517 DO i=1,n
518 x(i)=b(i)
519 END DO
520 DO i=1,n ! forward substitution
521 DO j=1,min(mp1-1,n-i)
522 x(j+i)=x(j+i)-w(j+1,i)*x(i)
523 END DO
524 END DO
525 DO i=n,1,-1 ! backward substitution
526 rxw=x(i)*w(1,i)
527 DO j=1,min(mp1-1,n-i)
528 rxw=rxw-w(j+1,i)*x(j+i)
529 END DO
530 x(i)=rxw
531 END DO
532 RETURN
533END SUBROUTINE dbcslv
534
541
542SUBROUTINE dbciel(w,mp1,n,v)
543 USE mpdef
544
545 implicit none
546 INTEGER(mpi) :: i
547 INTEGER(mpi) :: j
548 INTEGER(mpi) :: k
549 REAL(mpd) :: rxw
550
551 INTEGER(mpi), INTENT(IN) :: n
552 INTEGER(mpi), INTENT(IN) :: mp1
553 REAL(mpd), INTENT(IN) :: w(mp1,n)
554 REAL(mpd), INTENT(OUT) :: v(mp1,n)
555 ! N*M*(M-1) dot operations
556 DO i=n,1,-1
557 rxw=w(1,i)
558 DO j=i,max(1,i-mp1+1),-1
559 DO k=j+1,min(n,j+mp1-1)
560 rxw=rxw-v(1+abs(i-k),min(i,k))*w(1+k-j,j)
561 END DO
562 v(1+i-j,j)=rxw
563 rxw=0.0
564 END DO
565 END DO
566 RETURN
567END SUBROUTINE dbciel
568
575
576SUBROUTINE dbcinb(w,mp1,n, vs)
577 USE mpdef
578
579 implicit none
580 INTEGER(mpi) :: i
581 INTEGER(mpi) :: j
582 INTEGER(mpi) :: k
583 REAL(mpd) :: rxw
584
585 INTEGER(mpi), INTENT(IN) :: n
586 INTEGER(mpi), INTENT(IN) :: mp1
587 REAL(mpd), INTENT(IN) :: w(mp1,n)
588 REAL(mpd), INTENT(OUT) :: vs((n*n+n)/2)
589 ! N*M*(M-1) dot operations
590 DO i=n,1,-1
591 rxw=w(1,i)
592 DO j=i,max(1,i-mp1+1),-1
593 DO k=j+1,min(n,j+mp1-1)
594 rxw=rxw-vs((max(i,k)*(max(i,k)-1))/2+min(i,k))*w(1+k-j,j)
595 END DO
596 vs((i*i-i)/2+j)=rxw
597 rxw=0.0
598 END DO
599 ! DO J=MAX(1,I-MP1+1)-1,1,-1
600 ! VS((I*I-I)/2+J)=0.0
601 ! END DO
602 END DO
603 RETURN
604END SUBROUTINE dbcinb
605
612
613SUBROUTINE dbcinv(w,mp1,n, vs)
614 USE mpdef
615
616 implicit none
617 INTEGER(mpi) :: i
618 INTEGER(mpi) :: j
619 INTEGER(mpi) :: k
620 REAL(mpd) :: rxw
621
622 INTEGER(mpi), INTENT(IN) :: n
623 INTEGER(mpi), INTENT(IN) :: mp1
624 REAL(mpd), INTENT(IN) :: w(mp1,n)
625 REAL(mpd), INTENT(OUT) :: vs((n*n+n)/2)
626 ! N*N/2*(M-1) dot operations
627 DO i=n,1,-1
628 rxw=w(1,i)
629 DO j=i,1,-1
630 DO k=j+1,min(n,j+mp1-1)
631 rxw=rxw-vs((max(i,k)*(max(i,k)-1))/2+min(i,k))*w(1+k-j,j)
632 END DO
633 vs((i*i-i)/2+j)=rxw
634 rxw=0.0
635 END DO
636 END DO
637 RETURN
638END SUBROUTINE dbcinv
639
646
647SUBROUTINE dbcprv(w,mp1,n,b)
648 USE mpdef
649
650 implicit none
651 INTEGER(mpi) :: i
652 INTEGER(mpi) :: j
653 INTEGER(mpi) :: nj
654 REAL(mps) :: rho
655
656
657 INTEGER(mpi), INTENT(IN) :: n
658 INTEGER(mpi), INTENT(IN) :: mp1
659 REAL(mpd), INTENT(IN OUT) :: w(mp1,n)
660 REAL(mpd), INTENT(OUT) :: b(n)
661
662 REAL(mpd) :: ERR
663 INTEGER(mpi) :: irho(5)
664 ! ...
665 WRITE(*,*) ' '
666 WRITE(*,101)
667
668 DO i=1,n
669 err=sqrt(w(1,i))
670 nj=0
671 DO j=2,mp1
672 IF(i+1-j > 0.AND.nj < 5) THEN
673 nj=nj+1
674 rho=real(w(j,i+1-j)/(err*sqrt(w(1,i+1-j))),mps)
675 irho(nj)=nint(100.0*abs(rho),mpi)
676 IF(rho < 0.0) irho(nj)=-irho(nj)
677 END IF
678 END DO
679 WRITE(*,102) i,b(i),err,(irho(j),j=1,nj)
680 END DO
681 WRITE(*,103)
682101 FORMAT(5x,'i Param',7x,'error',7x,' c(i,i-1) c(i,i-2)'/)
683102 FORMAT(1x,i5,2g12.4,1x,5i9)
684103 FORMAT(33x,'(correlation coefficients in percent)')
685END SUBROUTINE dbcprv
686
692
693SUBROUTINE dbcprb(w,mp1,n)
694 USE mpdef
695
696 implicit none
697 INTEGER(mpi) :: i
698 INTEGER(mpi) :: j
699
700
701 INTEGER(mpi), INTENT(IN) :: mp1
702 INTEGER(mpi), INTENT(IN) :: n
703 REAL(mpd), INTENT(IN OUT) :: w(mp1,n)
704
705
706 ! ...
707 IF(mp1 > 6) RETURN
708 WRITE(*,*) ' '
709 WRITE(*,101)
710 DO i=1,n
711 WRITE(*,102) i,(w(j,i),j=1,mp1)
712 END DO
713 WRITE(*,*) ' '
714101 FORMAT(5x,'i Diag ')
715102 FORMAT(1x,i5,6g12.4)
716END SUBROUTINE dbcprb
717
718
719! (3) Symmetric band matrix of band width m=1: decomposition,
720! solution, complete, inverse and band part of the inverse
721
743SUBROUTINE db2dec(w,n, aux)
744 USE mpdef
745
746 implicit none
747 INTEGER(mpi) :: i
748 REAL(mpd) :: rxw
749
750 INTEGER(mpi), INTENT(IN OUT) :: n
751 REAL(mpd), INTENT(IN OUT) :: w(2,n)
752 REAL(mpd), INTENT(OUT) :: aux(n)
753
754 DO i=1,n
755 aux(i)=16.0_mpd*w(1,i) ! save diagonal elements
756 END DO
757 DO i=1,n-1
758 IF(w(1,i)+aux(i) /= aux(i)) THEN
759 w(1,i)=1.0_mpd/w(1,i)
760 rxw=w(2,i)*w(1,i)
761 w(1,i+1)=w(1,i+1)-w(2,i)*rxw
762 w(2,i)=rxw
763 ELSE ! singular
764 w(1,i)=0.0_mpd
765 w(2,i)=0.0_mpd
766 END IF
767 END DO
768 IF(w(1,n)+aux(n) > aux(n)) THEN ! N
769 w(1,n)=1.0_mpd/w(1,n)
770 ELSE ! singular
771 w(1,n)=0.0_mpd
772 END IF
773 RETURN
774END SUBROUTINE db2dec
775
782
783
784SUBROUTINE db2slv(w,n,b, x)
785 USE mpdef
786
787 implicit none
788 INTEGER(mpi) :: i
789
790 INTEGER(mpi), INTENT(IN) :: n
791 REAL(mpd), INTENT(IN) :: w(2,n)
792 REAL(mpd), INTENT(IN) :: b(n)
793 REAL(mpd), INTENT(OUT) :: x(n)
794
795 ! The equation W(original)*X=B is solved for X; input is B in vector X.
796 DO i=1,n
797 x(i)=b(i)
798 END DO
799 DO i=1,n-1 ! by forward substitution
800 x(i+1)=x(i+1)-w(2,i)*x(i)
801 END DO
802 x(n)=x(n)*w(1,n) ! and backward substitution
803 DO i=n-1,1,-1
804 x(i)=x(i)*w(1,i)-w(2,i)*x(i+1)
805 END DO
806 RETURN
807END SUBROUTINE db2slv
808
814
815SUBROUTINE db2iel(w,n, v)
816 USE mpdef
817
818 implicit none
819 INTEGER(mpi) :: i
820
821 INTEGER(mpi), INTENT(IN) :: n
822 REAL(mpd), INTENT(IN) :: w(2,n)
823 REAL(mpd), INTENT(OUT) :: v(2,n)
824
825 ! The band elements of the inverse of W(original) are calculated,
826 ! and stored in V in the same order as in W.
827 ! Remaining elements of the inverse are not calculated.
828 v(1,n )= w(1,n)
829 v(2,n-1)=-v(1,n )*w(2,n-1)
830 DO i=n-1,3,-1
831 v(1,i )= w(1,i )-v(2,i )*w(2,i )
832 v(2,i-1)=-v(1,i )*w(2,i-1)
833 END DO
834 v(1,2)= w(1,2)-v(2,2)*w(2,2)
835 v(2,1)=-v(1,2)*w(2,1)
836 v(1,1)= w(1,1)-v(2,1)*w(2,1)
837 RETURN
838END SUBROUTINE db2iel
839
840
841! (4) Symmetric band matrix of band width m=2: decomposition,
842! solution, complete, inverse and band part of the inverse
843
864
865SUBROUTINE db3dec(w,n, aux) ! decomposition (M=2)
866 USE mpdef
867
868 implicit none
869 INTEGER(mpi) :: i
870 REAL(mpd) :: rxw
871
872 INTEGER(mpi), INTENT(IN OUT) :: n
873 REAL(mpd), INTENT(IN OUT) :: w(3,n)
874 REAL(mpd), INTENT(OUT) :: aux(n)
875 ! decompos
876
877 DO i=1,n
878 aux(i)=16.0_mpd*w(1,i) ! save diagonal elements
879 END DO
880 DO i=1,n-2
881 IF(w(1,i)+aux(i) /= aux(i)) THEN
882 w(1,i)=1.0_mpd/w(1,i)
883 rxw=w(2,i)*w(1,i)
884 w(1,i+1)=w(1,i+1)-w(2,i)*rxw
885 w(2,i+1)=w(2,i+1)-w(3,i)*rxw
886 w(2,i)=rxw
887 rxw=w(3,i)*w(1,i)
888 w(1,i+2)=w(1,i+2)-w(3,i)*rxw
889 w(3,i)=rxw
890 ELSE ! singular
891 w(1,i)=0.0_mpd
892 w(2,i)=0.0_mpd
893 w(3,i)=0.0_mpd
894 END IF
895 END DO
896 IF(w(1,n-1)+aux(n-1) > aux(n-1)) THEN
897 w(1,n-1)=1.0_mpd/w(1,n-1)
898 rxw=w(2,n-1)*w(1,n-1)
899 w(1,n)=w(1,n)-w(2,n-1)*rxw
900 w(2,n-1)=rxw
901 ELSE ! singular
902 w(1,n-1)=0.0_mpd
903 w(2,n-1)=0.0_mpd
904 END IF
905 IF(w(1,n)+aux(n) > aux(n)) THEN
906 w(1,n)=1.0_mpd/w(1,n)
907 ELSE ! singular
908 w(1,n)=0.0_mpd
909 END IF
910 RETURN
911END SUBROUTINE db3dec
912
919
920SUBROUTINE db3slv(w,n,b, x)
921 USE mpdef
922
923 implicit none
924 INTEGER(mpi) :: i
925
926 INTEGER(mpi), INTENT(IN) :: n
927 REAL(mpd), INTENT(IN) :: w(3,n)
928 REAL(mpd), INTENT(IN) :: b(n)
929 REAL(mpd), INTENT(OUT) :: x(n)
930
931 DO i=1,n
932 x(i)=b(i)
933 END DO
934 DO i=1,n-2 ! by forward substitution
935 x(i+1)=x(i+1)-w(2,i)*x(i)
936 x(i+2)=x(i+2)-w(3,i)*x(i)
937 END DO
938 x(n)=x(n)-w(2,n-1)*x(n-1)
939 x(n)=x(n)*w(1,n) ! and backward substitution
940 x(n-1)=x(n-1)*w(1,n-1)-w(2,n-1)*x(n)
941 DO i=n-2,1,-1
942 x(i)=x(i)*w(1,i)-w(2,i)*x(i+1)-w(3,i)*x(i+2)
943 END DO
944 RETURN
945END SUBROUTINE db3slv
946
952
953SUBROUTINE db3iel(w,n, v)
954 USE mpdef
955
956 implicit none
957 INTEGER(mpi) :: i
958
959 INTEGER(mpi), INTENT(IN) :: n
960 REAL(mpd), INTENT(IN) :: w(3,n)
961 REAL(mpd), INTENT(OUT) :: v(3,n)
962
963 ! The band elements of the inverse of W(original) are calculated,
964 ! and stored in V in the same order as in W.
965 ! Remaining elements of the inverse are not calculated.
966 v(1,n )= w(1,n)
967 v(2,n-1)=-v(1,n )*w(2,n-1)
968 v(3,n-2)=-v(2,n-1)*w(2,n-2)-v(1,n )*w(3,n-2)
969 v(1,n-1)= w(1,n-1) -v(2,n-1)*w(2,n-1)
970 v(2,n-2)=-v(1,n-1)*w(2,n-2)-v(2,n-1)*w(3,n-2)
971 v(3,n-3)=-v(2,n-2)*w(2,n-3)-v(1,n-1)*w(3,n-3)
972 DO i=n-2,3,-1
973 v(1,i )= w(1,i ) -v(2,i )*w(2,i )-v(3,i)*w(3,i )
974 v(2,i-1)=-v(1,i )*w(2,i-1)-v(2,i)*w(3,i-1)
975 v(3,i-2)=-v(2,i-1)*w(2,i-2)-v(1,i)*w(3,i-2)
976 END DO
977 v(1,2)= w(1,2) -v(2,2)*w(2,2)-v(3,2)*w(3,2)
978 v(2,1)=-v(1,2)*w(2,1)-v(2,2)*w(3,1)
979 v(1,1)= w(1,1) -v(2,1)*w(2,1)-v(3,1)*w(3,1)
980 RETURN
981END SUBROUTINE db3iel
982
983
984! (5) Symmetric matrix with band structure, bordered by full row/col
985! - is not yet included -
986
987! SUBROUTINE BSOLV1(N,CU,RU,CK,RK,CH, BK,UH, AU) ! 1
988! Input: CU = 3*N array replaced by decomposition
989! RU N array rhs
990! CK diagonal element
991! RK rhs
992! CH N-vector
993
994! Aux: AU N-vector auxliliary array
995
996! Result: FK curvature
997! BK variance
998! UH smoothed data points
999
1000
1001! DOUBLE PRECISION CU(3,N),CI(3,N),CK,BK,AU(N),UH(N)
1002! ...
1003! CALL BDADEC(CU,3,N, AU) ! decomposition
1004! CALL DBASLV(CU,3,N, RU,UH) ! solve for zero curvature
1005! CALL DBASLV(CU,3,N, CH,AU) ! solve for aux. vector
1006! CTZ=0.0D0
1007! ZRU=0.0D0
1008! DO I=1,N
1009! CTZ=CTZ+CH(I)*AU(I) ! cT z
1010! ZRU=ZRU+RY(I)*AU(I) ! zT ru
1011! END DO
1012! BK=1.0D0/(CK-CTZ) ! variance of curvature
1013! FK=BK *(RK-ZRU) ! curvature
1014! DO I=1,N
1015! UH(I)=UH(I)-FK*AU(I) ! smoothed data points
1016! END DO
1017! RETURN
1018
1019! ENTRY BINV1(N,CU,CI, FK,AU)
1020! DOUBLE PRECISION CI(3,N)
1021! ...
1022! CALL DBAIBM(CU,3,N, CI) ! block part of inverse
1023! DO I=1,N
1024! CI(1,I)=CI(1,I)+FK*AU(I)*AU(I) ! diagonal elements
1025! IF(I.LT.N) CI(2,I)=CI(2,I)+FK*AU(I)*AU(I+1) ! next diagonal
1026! IF(I.LT.N-1) CI(3,I)=CI(3,I)+FK*AU(I)*AU(I+2) ! next diagonal
1027! END DO
1028
1029! END
1030
1039
1040SUBROUTINE dcfdec(w,n)
1041 USE mpdef
1042
1043 IMPLICIT NONE
1044 INTEGER(mpi), INTENT(IN) :: n
1045 REAL(mpd), INTENT(OUT) :: w((n*n+n)/2)
1046 INTEGER(mpi) :: i,j,k
1047 REAL(mpd) :: epsm,gamm,xchi,beta,delta,theta
1048
1049 epsm=epsilon(epsm) ! machine precision
1050 gamm=0.0_mpd ! max diagonal element
1051 xchi=0.0_mpd ! max off-diagonal element
1052 DO k=1,n
1053 gamm=max(gamm,abs(w((k*k+k)/2)))
1054 DO j=k+1,n
1055 xchi=max(xchi,abs(w((j*j-j)/2+k)))
1056 END DO
1057 END DO
1058 beta=sqrt(max(gamm,xchi/max(1.0_mpd,sqrt(real(n*n-1,mpd))),epsm))
1059 delta=epsm*max(1.0_mpd,gamm+xchi)
1060
1061 DO k=1,n
1062 DO i=1,k-1
1063 w((k*k-k)/2+i)=w((k*k-k)/2+i)*w((i*i+i)/2)
1064 END DO
1065 DO j=k+1,n
1066 DO i=1,k-1
1067 w((j*j-j)/2+k)=w((j*j-j)/2+k)-w((k*k-k)/2+i)*w((j*j-j)/2+i)
1068 END DO
1069 END DO
1070 theta=0.0_mpd
1071 DO j=k+1,n
1072 theta=max(theta,abs(w((j*j-j)/2+k)))
1073 END DO
1074 w((k*k+k)/2)=1.0_mpd/max(abs(w((k*k+k)/2)),(theta/beta)**2,delta)
1075 DO j=k+1,n
1076 w((j*j+j)/2)=w((j*j+j)/2)-w((j*j-j)/2+k)**2*w((k*k+k)/2)
1077 END DO
1078 END DO ! K
1079
1080END SUBROUTINE dcfdec
1081
1091
1092SUBROUTINE dbfdec(w,mp1,n)
1093 USE mpdef
1094
1095 IMPLICIT NONE
1096 INTEGER(mpi), INTENT(IN OUT) :: mp1
1097 INTEGER(mpi), INTENT(IN) :: n
1098 REAL(mpd), INTENT(OUT) :: w(mp1,n)
1099 INTEGER(mpi) :: i,j,k
1100 REAL(mpd) :: epsm,gamm,xchi,beta,delta,theta
1101
1102 epsm=epsilon(epsm) ! machine precision
1103 gamm=0.0_mpd ! max diagonal element
1104 xchi=0.0_mpd ! max off-diagonal element
1105 DO k=1,n
1106 gamm=max(gamm,abs(w(1,k)))
1107 DO j=2,min(mp1,n-k+1)
1108 xchi=max(xchi,abs(w(j,k)))
1109 END DO
1110 END DO
1111 beta=sqrt(max(gamm,xchi/max(1.0_mpd,sqrt(real(n*n-1,mpd))),epsm))
1112 delta=epsm*max(1.0_mpd,gamm+xchi)
1113
1114 DO k=1,n
1115 DO i=2,min(mp1,k)
1116 w(i,k-i+1)=w(i,k-i+1)*w(1,k-i+1)
1117 END DO
1118 DO j=2,min(mp1,n-k+1)
1119 DO i=max(2,j+k+1-mp1),k
1120 w(j,k)=w(j,k)-w(k-i+2,i-1)*w(j-i+k+1,i-1)
1121 END DO
1122 END DO
1123 theta=0.0_mpd
1124 DO j=2,min(mp1,n-k+1)
1125 theta=max(theta,abs(w(j,k)))
1126 END DO
1127 w(1,k)=1.0_mpd/max(abs(w(1,k)),(theta/beta)**2,delta)
1128 DO j=2,min(mp1,n-k+1)
1129 w(1,k+j-1)=w(1,k+j-1)-w(1,k)*w(j,k)**2
1130 END DO
1131 END DO ! K
1132
1133END SUBROUTINE dbfdec
1134
1135
dbcslv
subroutine dbcslv(w, mp1, n, b, x)
solution B -> X
Definition: Dbandmatrix.f90:504
dchslv
subroutine dchslv(w, n, b, x)
solution B -> X
Definition: Dbandmatrix.f90:296
dbcprv
subroutine dbcprv(w, mp1, n, b)
Print corr band matrix and pars.
Definition: Dbandmatrix.f90:648
dbcinv
subroutine dbcinv(w, mp1, n, vs)
VS = inverse symmetric matrix.
Definition: Dbandmatrix.f90:614
dbfdec
subroutine dbfdec(w, mp1, n)
Decomposition of symmetric band matrix.
Definition: Dbandmatrix.f90:1093
dbcinb
subroutine dbcinb(w, mp1, n, vs)
VS = band part of inverse symmetric matrix.
Definition: Dbandmatrix.f90:577
condes
real(mps) function condes(w, n, aux)
Etimate condition.
Definition: Dbandmatrix.f90:385
db3iel
subroutine db3iel(w, n, v)
V = inverse band matrix elements.
Definition: Dbandmatrix.f90:954
dbcdec
subroutine dbcdec(w, mp1, n, aux)
Decomposition of symmetric band matrix.
Definition: Dbandmatrix.f90:460
dcfdec
subroutine dcfdec(w, n)
Decomposition of symmetric matrix.
Definition: Dbandmatrix.f90:1041
db2dec
subroutine db2dec(w, n, aux)
Decomposition (M=1).
Definition: Dbandmatrix.f90:744
db3dec
subroutine db3dec(w, n, aux)
Decomposition (M=2).
Definition: Dbandmatrix.f90:866
dbcprb
subroutine dbcprb(w, mp1, n)
Print band matrix.
Definition: Dbandmatrix.f90:694
dbciel
subroutine dbciel(w, mp1, n, v)
V = inverse band matrix elements.
Definition: Dbandmatrix.f90:543
db2slv
subroutine db2slv(w, n, b, x)
solution B -> X
Definition: Dbandmatrix.f90:785
db3slv
subroutine db3slv(w, n, b, x)
solution B -> X
Definition: Dbandmatrix.f90:921
dchdec
subroutine dchdec(w, n, aux)
Decomposition of symmetric matrix.
Definition: Dbandmatrix.f90:246
dchinv
subroutine dchinv(w, n, v)
inversion
Definition: Dbandmatrix.f90:346
db2iel
subroutine db2iel(w, n, v)
V = inverse band matrix elements.
Definition: Dbandmatrix.f90:816
mpdef
Definition of constants.
Definition: mpdef.f90:24
mpdef::mpd
integer, parameter mpd
double precision
Definition: mpdef.f90:38
mpdef::mps
integer, parameter mps
single precision
Definition: mpdef.f90:37
mpdef::mpi
integer, parameter mpi
integer
Definition: mpdef.f90:35
Generated on Thu Dec 5 2024 14:18:59 for Millepede-II by doxygen 1.9.3