SND@LHC Software
Loading...
Searching...
No Matches
mpnum.f90
Go to the documentation of this file.
1
2! Code converted using TO_F90 by Alan Miller
3! Date: 2012-03-16 Time: 11:07:48
4
72
73!----------------------------------------------------------------------
92
93SUBROUTINE sqminv(v,b,n,nrank,diag,next) ! matrix inversion
94 USE mpdef
95
96 IMPLICIT NONE
97 INTEGER(mpi) :: i
98 INTEGER(mpi) :: ij
99 INTEGER(mpi) :: j
100 INTEGER(mpi) :: jj
101 INTEGER(mpi) :: jk
102 INTEGER(mpi) :: jl
103 INTEGER(mpi) :: k
104 INTEGER(mpi) :: kk
105 INTEGER(mpi) :: l
106 INTEGER(mpi) :: last
107 INTEGER(mpi) :: lk
108 INTEGER(mpi) :: next0
109
110 REAL(mpd), INTENT(IN OUT) :: v(*)
111 REAL(mpd), INTENT(OUT) :: b(n)
112 INTEGER(mpi), INTENT(IN) :: n
113 INTEGER(mpi), INTENT(OUT) :: nrank
114 REAL(mpd), INTENT(OUT) :: diag(n)
115 INTEGER(mpi), INTENT(OUT) :: next(n)
116 REAL(mpd) :: vkk
117 REAL(mpd) :: vjk
118
119 !REAL(mpd), PARAMETER :: eps=1.0E-10_mpd
120 REAL(mpd) eps
121 ! ...
122 eps = 16.0_mpd * epsilon(eps) ! 16 * precision(mpd)
123
124 next0=1
125 l=1
126 DO i=1,n
127 next(i)=i+1 ! set "next" pointer
128 diag(i)=abs(v((i*i+i)/2)) ! save abs of diagonal elements
129 END DO
130 next(n)=-1 ! end flag
131
132 nrank=0
133 DO i=1,n ! start of loop
134 k =0
135 vkk=0.0_mpd
136
137 j=next0
138 last=0
13905 IF(j > 0) THEN
140 jj=(j*j+j)/2
141 IF(abs(v(jj)) > max(abs(vkk),eps*diag(j))) THEN
142 vkk=v(jj)
143 k=j
144 l=last
145 END IF
146 last=j
147 j=next(last)
148 GO TO 05
149 END IF
150
151 IF(k /= 0) THEN ! pivot found
152 kk=(k*k+k)/2
153 IF(l == 0) THEN
154 next0=next(k)
155 ELSE
156 next(l)=next(k)
157 END IF
158 next(k)=0 ! index is used, reset
159 nrank=nrank+1 ! increase rank and ...
160 vkk =1.0/vkk
161 v(kk) =-vkk
162 b(k) =b(k)*vkk
163 jk =kk-k
164 jl =0
165 DO j=1,n ! elimination
166 IF(j == k) THEN
167 jk=kk
168 jl=jl+j
169 ELSE
170 IF(j < k) THEN
171 jk=jk+1
172 ELSE
173 jk=jk+j-1
174 END IF
175 vjk =v(jk)
176 v(jk)=vkk*vjk
177 b(j) =b(j)-b(k)*vjk
178 lk =kk-k
179 DO l=1,j
180 jl=jl+1
181 IF(l == k) THEN
182 lk=kk
183 ELSE
184 IF(l < k) THEN
185 lk=lk+1
186 ELSE
187 lk=lk+l-1
188 END IF
189 v(jl)=v(jl)-v(lk)*vjk
190 END IF
191 END DO
192 END IF
193 END DO
194 ELSE
195 DO k=1,n
196 IF(next(k) /= 0) THEN
197 b(k)=0.0_mpd ! clear vector element
198 DO j=1,k
199 IF(next(j) /= 0) v((k*k-k)/2+j)=0.0_mpd ! clear matrix row/col
200 END DO
201 END IF
202 END DO
203 GO TO 10
204 END IF
205 END DO ! end of loop
206 10 DO ij=1,(n*n+n)/2
207 v(ij)=-v(ij) ! finally reverse sign of all matrix elements
208 END DO
209END SUBROUTINE sqminv
210
223
224SUBROUTINE sqminl(v,b,n,nrank,diag,next) !
225 USE mpdef
226
227 IMPLICIT NONE
228 INTEGER(mpi) :: i
229 INTEGER(mpi) :: j
230 INTEGER(mpi) :: k
231 INTEGER(mpi) :: l
232 INTEGER(mpi) :: last
233 INTEGER(mpi) :: next0
234
235 REAL(mpd), INTENT(IN OUT) :: v(*)
236 REAL(mpd), INTENT(OUT) :: b(n)
237 INTEGER(mpi), INTENT(IN) :: n
238 INTEGER(mpi), INTENT(OUT) :: nrank
239 REAL(mpd), INTENT(OUT) :: diag(n)
240 INTEGER(mpi), INTENT(OUT) :: next(n)
241 INTEGER(mpl) :: i8
242 INTEGER(mpl) :: j8
243 INTEGER(mpl) :: jj
244 INTEGER(mpl) :: k8
245 INTEGER(mpl) :: kk
246 INTEGER(mpl) :: kkmk
247 INTEGER(mpl) :: jk
248 INTEGER(mpl) :: jl
249 INTEGER(mpl) :: llk
250 INTEGER(mpl) :: ljl
251
252 REAL(mpd) :: vkk
253 REAL(mpd) :: vjk
254
255 REAL(mpd), PARAMETER :: eps=1.0e-10_mpd
256 ! ...
257 next0=1
258 l=1
259 DO i=1,n
260 i8=int8(i)
261 next(i)=i+1 ! set "next" pointer
262 diag(i)=abs(v((i8*i8+i8)/2)) ! save abs of diagonal elements
263 END DO
264 next(n)=-1 ! end flag
265
266 nrank=0
267 DO i=1,n ! start of loop
268 k =0
269 vkk=0.0_mpd
270 j=next0
271 last=0
27205 IF(j > 0) THEN
273 j8=int8(j)
274 jj=(j8*j8+j8)/2
275 IF(abs(v(jj)) > max(abs(vkk),eps*diag(j))) THEN
276 vkk=v(jj)
277 k=j
278 l=last
279 END IF
280 last=j
281 j=next(last)
282 GO TO 05
283 END IF
284
285 IF(k /= 0) THEN ! pivot found
286 k8=int8(k)
287 kk=(k8*k8+k8)/2
288 kkmk=kk-k8
289 IF(l == 0) THEN
290 next0=next(k)
291 ELSE
292 next(l)=next(k)
293 END IF
294 next(k)=0 ! index is used, reset
295 nrank=nrank+1 ! increase rank and ...
296 vkk =1.0/vkk
297 v(kk) =-vkk
298 b(k) =b(k)*vkk
299 ! elimination
300 jk =kkmk
301 DO j=1,n
302 IF(j == k) THEN
303 jk=kk
304 ELSE
305 IF(j < k) THEN
306 jk=jk+1
307 ELSE
308 jk=jk+int8(j)-1
309 END IF
310 v(jk)=v(jk)*vkk
311 END IF
312 END DO
313 ! parallelize row loop
314 ! slot of 128 'J' for next idle thread
315 !$OMP PARALLEL DO &
316 !$OMP PRIVATE(JL,JK,L,LJL,LLK,VJK,J8) &
317 !$OMP SCHEDULE(DYNAMIC,128)
318 DO j=n,1,-1
319 j8=int8(j)
320 jl=j8*(j8-1)/2
321 IF(j /= k) THEN
322 IF(j < k) THEN
323 jk=kkmk+j8
324 ELSE
325 jk=k8+jl
326 END IF
327 vjk =v(jk)/vkk
328 b(j) =b(j)-b(k)*vjk
329 ljl=jl
330 llk=kkmk
331 DO l=1,min(j,k-1)
332 ljl=ljl+1
333 llk=llk+1
334 v(ljl)=v(ljl)-v(llk)*vjk
335 END DO
336 ljl=ljl+1
337 llk=kk
338 DO l=k+1,j
339 ljl=ljl+1
340 llk=llk+l-1
341 v(ljl)=v(ljl)-v(llk)*vjk
342 END DO
343 END IF
344 END DO
345 !$OMP END PARALLEL DO
346 ELSE
347 DO k=1,n
348 k8=int8(k)
349 kk=(k8*k8-k8)/2
350 IF(next(k) /= 0) THEN
351 b(k)=0.0_mpd ! clear vector element
352 DO j=1,k
353 IF(next(j) /= 0) v(kk+int8(j))=0.0_mpd ! clear matrix row/col
354 END DO
355 END IF
356 END DO
357 GO TO 10
358 END IF
359 END DO ! end of loop
360 10 DO jj=1,(int8(n)*int8(n)+int8(n))/2
361 v(jj)=-v(jj) ! finally reverse sign of all matrix elements
362 END DO
363END SUBROUTINE sqminl
364
376
377SUBROUTINE devrot(n,diag,u,v,work,iwork) ! diagonalization
378 USE mpdef
379
380 IMPLICIT NONE
381
382 INTEGER(mpi), INTENT(IN) :: n
383 REAL(mpd), INTENT(OUT) :: diag(n)
384 REAL(mpd), INTENT(OUT) :: u(n,n)
385 REAL(mpd), INTENT(IN) :: v(*)
386 REAL(mpd), INTENT(OUT) :: work(n)
387 INTEGER(mpi), INTENT(OUT) :: iwork(n)
388
389
390 INTEGER(mpi), PARAMETER :: itmax=30
391 REAL(mpd), PARAMETER :: tol=epsilon(tol)
392 REAL(mpd), PARAMETER :: eps=epsilon(eps)
393
394 REAL(mpd) :: f
395 REAL(mpd) :: g
396 REAL(mpd) :: h
397 REAL(mpd) :: sh
398 REAL(mpd) :: hh
399 REAL(mpd) :: b
400 REAL(mpd) :: p
401 REAL(mpd) :: r
402 REAL(mpd) :: s
403 REAL(mpd) :: c
404 REAL(mpd) :: workd
405
406 INTEGER(mpi) :: ij
407 INTEGER(mpi) :: i
408 INTEGER(mpi) :: j
409 INTEGER(mpi) :: k
410 INTEGER(mpi) :: l
411 INTEGER(mpi) :: m
412 INTEGER(mpi) :: ll
413 ! ...
414 ! 1. part: symmetric matrix V reduced to tridiagonal from
415 ij=0
416 DO i=1,n
417 DO j=1,i
418 ij=ij+1
419 u(i,j)=v(ij) ! copy half of symmetric matirx
420 END DO
421 END DO
422
423 DO i=n,2,-1
424 l=i-2
425 f=u(i,i-1)
426 g=0.0_mpd
427 IF(l /= 0) THEN
428 DO k=1,l
429 IF(abs(u(i,k)) > tol) g=g+u(i,k)*u(i,k)
430 END DO
431 h=g+f*f
432 END IF
433 IF(g < tol) THEN ! G too small
434 work(i)=f ! skip transformation
435 h =0.0_mpd
436 ELSE
437 l=l+1
438 sh=sqrt(h)
439 IF(f >= 0.0_mpd) sh=-sh
440 g=sh
441 work(i)=sh
442 h=h-f*g
443 u(i,i-1)=f-g
444 f=0.0_mpd
445 DO j=1,l
446 u(j,i)=u(i,j)/h
447 g=0.0_mpd
448 ! form element of a u
449 DO k=1,j
450 IF(abs(u(j,k)) > tol.AND.abs(u(i,k)) > tol) THEN
451 g=g+u(j,k)*u(i,k)
452 END IF
453 END DO
454 DO k=j+1,l
455 IF(abs(u(k,j)) > tol.AND.abs(u(i,k)) > tol) THEN
456 g=g+u(k,j)*u(i,k)
457 END IF
458 END DO
459 work(j)=g/h
460 f=f+g*u(j,i)
461 END DO
462 ! form k
463 hh=f/(h+h)
464 ! form reduced a
465 DO j=1,l
466 f=u(i,j)
467 work(j)=work(j)-hh*f
468 g=work(j)
469 DO k=1,j
470 u(j,k)=u(j,k)-f*work(k)-g*u(i,k)
471 END DO
472 END DO
473 END IF
474 diag(i)=h
475 END DO
476
477 diag(1)=0.0_mpd
478 work(1)=0.0_mpd
479
480 ! accumulation of transformation matrices
481 DO i=1,n
482 IF(diag(i) /= 0.0) THEN
483 DO j=1,i-1
484 g=0.0_mpd
485 DO k=1,i-1
486 g=g+u(i,k)*u(k,j)
487 END DO
488 DO k=1,i-1
489 u(k,j)=u(k,j)-g*u(k,i)
490 END DO
491 END DO
492 END IF
493 diag(i)=u(i,i)
494 u(i,i)=1.0_mpd
495 DO j=1,i-1
496 u(i,j)=0.0_mpd
497 u(j,i)=0.0_mpd
498 END DO
499 END DO
500
501 ! 2. part: diagonalization of tridiagonal matrix
502 DO i=2,n
503 work(i-1)=work(i)
504 END DO
505 work(n)=0.0_mpd
506 b=0.0_mpd
507 f=0.0_mpd
508
509 DO l=1,n
510 j=0
511 h=eps*(abs(diag(l))+abs(work(l)))
512 IF(b < h) b=h
513 DO m=l,n
514 IF(abs(work(m)) <= b) GO TO 10 ! look for small sub-diagonal element
515 END DO
516 m=l
51710 IF(m == l) GO TO 30
518 ! next iteration
51920 IF(j == itmax) THEN
520 WRITE(*,*) 'DEVROT: Iteration limit reached'
521 CALL peend(32,'Aborted, iteration limit reached in diagonalization')
522 stop
523 END IF
524 j=j+1
525 g=diag(l)
526 p=(diag(l+1)-g)/(2.0_mpd*work(l))
527 r=sqrt(1.0_mpd+p*p)
528 diag(l)=work(l)
529 IF(p < 0.0_mpd) diag(l)=diag(l)/(p-r)
530 IF(p >= 0.0_mpd) diag(l)=diag(l)/(p+r)
531 h=g-diag(l)
532 DO i=l+1,n
533 diag(i)=diag(i)-h
534 END DO
535 f=f+h
536 ! QL transformation
537 p=diag(m)
538 c=1.0_mpd
539 s=0.0_mpd
540 DO i=m-1,l,-1 ! reverse loop
541 g=c*work(i)
542 h=c*p
543 IF(abs(p) >= abs(work(i))) THEN
544 c=work(i)/p
545 r=sqrt(1.0_mpd+c*c)
546 work(i+1)=s*p*r
547 s=c/r
548 c=1.0_mpd/r
549 ELSE
550 c=p/work(i)
551 r=sqrt(1.0_mpd+c*c)
552 work(i+1)=s*work(i)*r
553 s=1.0_mpd/r
554 c=c/r
555 END IF
556 p=c*diag(i)-s*g
557 diag(i+1)=h+s*(c*g+s*diag(i))
558 ! form vector
559 DO k=1,n
560 h=u(k,i+1)
561 u(k,i+1)=s*u(k,i)+c*h
562 u(k,i)=c*u(k,i)-s*h
563 END DO
564 END DO
565 work(l)=s*p
566 diag(l)=c*p
567 IF(abs(work(l)) > b) GO TO 20 ! next iteration
56830 diag(l)=diag(l)+f
569 END DO
570 DO i=1,n
571 iwork(i)=i
572 END DO
573
574 m=1
57540 m=1+3*m ! determine initial increment
576 IF(m <= n) GO TO 40
57750 m=m/3
578 DO j=1,n-m ! sort with increment M
579 l=j
58060 IF(diag(iwork(l+m)) > diag(iwork(l))) THEN ! compare
581 ll=iwork(l+m) ! exchange the two index values
582 iwork(l+m)=iwork(l)
583 iwork(l)=ll
584 l=l-m
585 IF(l > 0) GO TO 60
586 END IF
587 END DO
588 IF(m > 1) GO TO 50
589
590 DO i=1,n
591 IF(iwork(i) /= i) THEN
592 ! move vector from position I to the work area
593 workd=diag(i)
594 DO l=1,n
595 work(l)=u(l,i)
596 END DO
597 k=i
59870 j=k
599 k=iwork(j)
600 iwork(j)=j
601 IF(k /= i) THEN
602 ! move vector from position K to the (free) position J
603 diag(j)=diag(k)
604 DO l=1,n
605 u(l,j)=u(l,k)
606 END DO
607 GO TO 70
608 END IF
609 ! move vector from the work area to position J
610 diag(j)=workd
611 DO l=1,n
612 u(l,j)=work(l)
613 END DO
614 END IF
615 END DO
616END SUBROUTINE devrot
617
619SUBROUTINE devsig(n,diag,u,b,coef)
620 USE mpdef
621
622 IMPLICIT NONE
623
624 INTEGER(mpi), INTENT(IN) :: n
625 REAL(mpd), INTENT(IN) :: diag(n)
626 REAL(mpd), INTENT(IN) :: u(n,n)
627 REAL(mpd), INTENT(IN) :: b(n)
628 REAL(mpd), INTENT(OUT) :: coef(n)
629 INTEGER(mpi) :: i
630 INTEGER(mpi) :: j
631 REAL(mpd) :: dsum
632 ! ...
633 DO i=1,n
634 coef(i)=0.0_mpd
635 IF(diag(i) > 0.0_mpd) THEN
636 dsum=0.0_mpd
637 DO j=1,n
638 dsum=dsum+u(j,i)*b(j)
639 END DO
640 coef(i)=abs(dsum)/sqrt(diag(i))
641 END IF
642 END DO
643END SUBROUTINE devsig
644
645
656
657SUBROUTINE devsol(n,diag,u,b,x,work)
658 USE mpdef
659
660 IMPLICIT NONE
661
662 INTEGER(mpi), INTENT(IN) :: n
663 REAL(mpd), INTENT(IN) :: diag(n)
664 REAL(mpd), INTENT(IN) :: u(n,n)
665 REAL(mpd), INTENT(IN) :: b(n)
666 REAL(mpd), INTENT(OUT) :: x(n)
667 REAL(mpd), INTENT(OUT) :: work(n)
668 INTEGER(mpi) :: i
669 INTEGER(mpi) :: j
670 INTEGER(mpi) :: jj
671 REAL(mpd) :: s
672 ! ...
673 DO j=1,n
674 s=0.0_mpd
675 work(j)=0.0_mpd
676 IF(diag(j) /= 0.0_mpd) THEN
677 DO i=1,n
678 ! j-th eigenvector is U(.,J)
679 s=s+u(i,j)*b(i)
680 END DO
681 work(j)=s/diag(j)
682 END IF
683 END DO
684
685 DO j=1,n
686 s=0.0_mpd
687 DO jj=1,n
688 s=s+u(j,jj)*work(jj)
689 END DO
690 x(j)=s
691 END DO
692! WRITE(*,*) 'DEVSOL'
693! WRITE(*,*) 'X ',X
694END SUBROUTINE devsol
695
703
704SUBROUTINE devinv(n,diag,u,v)
705 USE mpdef
706
707 IMPLICIT NONE
708 INTEGER(mpi) :: i
709 INTEGER(mpi) :: ij
710 INTEGER(mpi) :: j
711 INTEGER(mpi) :: k
712
713 INTEGER(mpi), INTENT(IN) :: n
714 REAL(mpd), INTENT(IN) :: diag(n)
715 REAL(mpd), INTENT(IN) :: u(n,n)
716 REAL(mpd), INTENT(OUT) :: v(*)
717 REAL(mpd) :: dsum
718 ! ...
719 ij=0
720 DO i=1,n
721 DO j=1,i
722 ij=ij+1
723 dsum=0.0_mpd
724 DO k=1,n
725 IF(diag(k) /= 0.0_mpd) THEN
726 dsum=dsum+u(i,k)*u(j,k)/diag(k)
727 END IF
728 END DO
729 v(ij)=dsum
730 END DO
731 END DO
732END SUBROUTINE devinv
733
734
753SUBROUTINE choldc(g,n)
754 USE mpdef
755
756 IMPLICIT NONE
757 INTEGER(mpi) :: i
758 INTEGER(mpi) :: ii
759 INTEGER(mpi) :: j
760 INTEGER(mpi) :: jj
761 INTEGER(mpi) :: k
762 INTEGER(mpi) :: kk
763
764 REAL(mpd), INTENT(IN OUT) :: g(*)
765 INTEGER(mpi), INTENT(IN) :: n
766
767 REAL(mpd) :: ratio
768 ! ...
769 ii=0
770 DO i=1,n
771 ii=ii+i
772 IF(g(ii) /= 0.0) g(ii)=1.0/g(ii) ! (I,I) div !
773 jj=ii
774 DO j=i+1,n
775 ratio=g(i+jj)*g(ii) ! (I,J) (I,I)
776 kk=jj
777 DO k=j,n
778 g(kk+j)=g(kk+j)-g(kk+i)*ratio ! (K,J) (K,I)
779 kk=kk+k
780 END DO ! K
781 g(i+jj)=ratio ! (I,J)
782 jj=jj+j
783 END DO ! J
784 END DO ! I
785 RETURN
786END SUBROUTINE choldc
787
799SUBROUTINE cholsl(g,x,n)
800 USE mpdef
801
802 IMPLICIT NONE
803 REAL(mpd) :: dsum
804 INTEGER(mpi) :: i
805 INTEGER(mpi) :: ii
806 INTEGER(mpi) :: k
807 INTEGER(mpi) :: kk
808
809 REAL(mpd), INTENT(IN) :: g(*)
810 REAL(mpd), INTENT(IN OUT) :: x(n)
811 INTEGER(mpi), INTENT(IN) :: n
812
813 ii=0
814 DO i=1,n
815 dsum=x(i)
816 DO k=1,i-1
817 dsum=dsum-g(k+ii)*x(k) ! (K,I)
818 END DO
819 x(i)=dsum
820 ii=ii+i
821 END DO
822 DO i=n,1,-1
823 dsum=x(i)*g(ii) ! (I,I)
824 kk=ii
825 DO k=i+1,n
826 dsum=dsum-g(kk+i)*x(k) ! (K,I)
827 kk=kk+k
828 END DO
829 x(i)=dsum
830 ii=ii-i
831 END DO
832 RETURN
833END SUBROUTINE cholsl
834
845SUBROUTINE cholin(g,v,n)
846 USE mpdef
847
848 IMPLICIT NONE
849 REAL(mpd) :: dsum
850 INTEGER(mpi) :: i
851 INTEGER(mpi) :: ii
852 INTEGER(mpi) :: j
853 INTEGER(mpi) :: k
854 INTEGER(mpi) :: l
855 INTEGER(mpi) :: m
856
857 REAL(mpd), INTENT(IN) :: g(*)
858 REAL(mpd), INTENT( OUT) :: v(*)
859 INTEGER(mpi), INTENT(IN) :: n
860
861 ii=(n*n-n)/2
862 DO i=n,1,-1
863 dsum=g(ii+i) ! (I,I)
864 DO j=i,1,-1
865 DO k=j+1,n
866 l=min(i,k)
867 m=max(i,k)
868 dsum=dsum-g(j+(k*k-k)/2)*v(l+(m*m-m)/2) ! (J,K) (I,K)
869 END DO
870 v(ii+j)=dsum ! (I,J)
871 dsum=0.0_mpd
872 END DO
873 ii=ii-i+1
874 END DO
875END SUBROUTINE cholin
876
877! variable band matrix operations ----------------------------------
878
899
900SUBROUTINE vabdec(n,val,ilptr)
901 USE mpdef
902
903 IMPLICIT NONE
904
905 INTEGER(mpi) :: i
906 INTEGER(mpi) :: in
907 INTEGER(mpi) :: j
908 INTEGER(mpi) :: k
909 INTEGER(mpi) :: kj
910 INTEGER(mpi) :: mj
911 INTEGER(mpi) :: mk
912 REAL(mpd) :: sn
913 REAL(mpd) :: beta
914 REAL(mpd) :: delta
915 REAL(mpd) :: theta
916
917 INTEGER(mpi), INTENT(IN) :: n
918 REAL(mpd), INTENT(IN OUT) :: val(*)
919 INTEGER(mpi), INTENT(IN) :: ilptr(n)
920
921 REAL(mpd) :: dgamma
922 REAL(mpd) :: xi
923 REAL(mpd) :: valkj
924
925 REAL(mpd), PARAMETER :: one=1.0_mpd
926 REAL(mpd), PARAMETER :: two=2.0_mpd
927 REAL(mpd), PARAMETER :: eps = epsilon(eps)
928
929 WRITE(*,*) 'Variable band matrix Cholesky decomposition'
930
931 dgamma=0.0_mpd
932 i=1
933 DO j=1,ilptr(n) ! loop thrugh all matrix elements
934 IF(ilptr(i) == j) THEN ! diagonal element
935 IF(val(j) <= 0.0_mpd) GO TO 01 ! exit loop for negative diag
936 dgamma=max(dgamma,abs(val(j))) ! max diagonal element
937 i=i+1
938 END IF
939 END DO
940 i=n+1
94101 in=i-1 ! IN positive diagonal elements
942 WRITE(*,*) ' ',in,' positive diagonal elements'
943 xi=0.0_mpd
944 i=1
945 DO j=1,ilptr(in) ! loop for positive diagonal elements
946 ! through all matrix elements
947 IF(ilptr(i) == j) THEN ! diagonal element
948 i=i+1
949 ELSE
950 xi=max(xi,abs(val(j))) ! Xi = abs(max) off-diagonal element
951 END IF
952 END DO
953
954 delta=eps*max(1.0_mpd,dgamma+xi)
955 sn=1.0_mpd
956 IF(n > 1) sn=1.0_mpd/sqrt(real(n*n-1,mpd))
957 beta=sqrt(max(eps,dgamma,xi*sn)) ! beta
958 WRITE(*,*) ' DELTA and BETA ',delta,beta
959
960 DO k=2,n
961 mk=k-ilptr(k)+ilptr(k-1)+1
962
963 theta=0.0_mpd
964
965 DO j=mk,k
966 mj=j-ilptr(j)+ilptr(j-1)+1
967 kj=ilptr(k)-k+j ! index kj
968
969 DO i=max(mj,mk),j-1
970 val(kj)=val(kj) & ! L_kj := L_kj - L_ki D_ii L_ji
971 -val(ilptr(k)-k+i)*val(ilptr(i))*val(ilptr(j)-j+i)
972
973 END DO !
974
975 theta=max(theta,abs(val(kj))) ! maximum value of row
976
977 IF(j /= k) THEN
978 IF(val(ilptr(j)) /= 0.0_mpd) THEN
979 val(kj)=val(kj)/val(ilptr(j))
980 ELSE
981 val(kj)=0.0_mpd
982 END IF
983 END IF ! L_kj := L_kj/D_jj ! D_kk
984
985 IF(j == k) THEN
986 valkj=val(kj)
987 IF(k <= in) THEN
988 val(kj)=max(abs(val(kj)),(theta/beta)**2,delta)
989 IF(valkj /= val(kj)) THEN
990 WRITE(*,*) ' Index K=',k
991 WRITE(*,*) ' ',valkj,val(kj), (theta/beta)**2,delta,theta
992 END IF
993 END IF
994 END IF
995 END DO ! J
996
997 END DO ! K
998
999 DO k=1,n
1000 IF(val(ilptr(k)) /= 0.0_mpd) val(ilptr(k))=1.0_mpd/val(ilptr(k))
1001 END DO
1002
1003 RETURN
1004END SUBROUTINE vabdec
1005
1011
1012SUBROUTINE vabmmm(n,val,ilptr)
1013 USE mpdef
1014
1015 IMPLICIT NONE
1016 INTEGER(mpi) :: k
1017 INTEGER(mpi) :: kp
1018 INTEGER(mpi) :: kr
1019 INTEGER(mpi) :: ks
1020 INTEGER(mpi), INTENT(IN) :: n
1021 REAL(mpd), INTENT(IN OUT) :: val(*)
1022 INTEGER(mpi), INTENT(IN) :: ilptr(n)
1023 kr=1
1024 ks=1
1025 kp=1
1026 DO k=1,n
1027 IF(val(ilptr(k)) > val(ilptr(ks))) ks=k
1028 IF(val(ilptr(k)) < val(ilptr(kr))) kr=k
1029 IF(val(ilptr(k)) > 0.0.AND.val(ilptr(k)) < val(ilptr(kp))) kp=k
1030 END DO
1031 WRITE(*,*) ' Index value ',ks,val(ilptr(ks))
1032 WRITE(*,*) ' Index value ',kp,val(ilptr(kp))
1033 WRITE(*,*) ' Index value ',kr,val(ilptr(kr))
1034
1035 RETURN
1036END SUBROUTINE vabmmm
1037
1048
1049SUBROUTINE vabslv(n,val,ilptr,x)
1050 USE mpdef
1051
1052 IMPLICIT NONE
1053 INTEGER(mpi) :: j
1054 INTEGER(mpi) :: k
1055 INTEGER(mpi) :: mk
1056
1057 INTEGER(mpi), INTENT(IN) :: n
1058 REAL(mpd), INTENT(IN OUT) :: val(*)
1059 INTEGER(mpi), INTENT(IN) :: ilptr(n)
1060 REAL(mpd), INTENT(IN OUT) :: x(n)
1061 ! ...
1062 DO k=1,n ! forward loop
1063 mk=k-ilptr(k)+ilptr(k-1)+1
1064 DO j=mk,k-1
1065 x(k)=x(k)-val(ilptr(k)-k+j)*x(j) ! X_k := X_k - L_kj B_j
1066 END DO
1067 END DO ! K
1068
1069 DO k=1,n ! divide by diagonal elements
1070 x(k)=x(k)*val(ilptr(k)) ! X_k := X_k*D_kk
1071 END DO
1072
1073 DO k=n,1,-1 ! backward loop
1074 mk=k-ilptr(k)+ilptr(k-1)+1
1075 DO j=mk,k-1
1076 x(j)=x(j)-val(ilptr(k)-k+j)*x(k) ! X_j := X_j - L_kj X_k
1077 END DO
1078 END DO ! K
1079END SUBROUTINE vabslv
1080
1081! matrix/vector products -------------------------------------------
1082
1089
1090REAL(mpd) FUNCTION dbdot(n,dx,dy)
1091 USE mpdef
1092
1093 IMPLICIT NONE
1094 INTEGER(mpi) :: i
1095 REAL(mpd) :: dtemp
1096
1097 INTEGER(mpi), INTENT(IN) :: n
1098 REAL(mpd), INTENT(IN) :: dx(*)
1099 REAL(mpd), INTENT(IN) :: dy(*)
1100 ! ...
1101 dtemp=0.0_mpd
1102 DO i = 1,mod(n,5)
1103 dtemp=dtemp+dx(i)*dy(i)
1104 END DO
1105 DO i =mod(n,5)+1,n,5
1106 dtemp=dtemp+dx(i)*dy(i)+dx(i+1)*dy(i+1)+dx(i+2)*dy(i+2) &
1107 +dx(i+3)*dy(i+3)+dx(i+4)*dy(i+4)
1108 END DO
1109 dbdot=dtemp
1110END FUNCTION dbdot
1111
1120
1121SUBROUTINE dbaxpy(n,da,dx,dy)
1122 USE mpdef
1123
1124 IMPLICIT NONE
1125 INTEGER(mpi) :: i
1126
1127 INTEGER(mpi), INTENT(IN) :: n
1128 REAL(mpd), INTENT(IN) :: dx(*)
1129 REAL(mpd), INTENT(IN OUT) :: dy(*)
1130 REAL(mpd), INTENT(IN) :: da
1131 ! ...
1132 DO i=1,mod(n,4)
1133 dy(i)=dy(i)+da*dx(i)
1134 END DO
1135 DO i=mod(n,4)+1,n,4
1136 dy(i )=dy(i )+da*dx(i )
1137 dy(i+1)=dy(i+1)+da*dx(i+1)
1138 dy(i+2)=dy(i+2)+da*dx(i+2)
1139 dy(i+3)=dy(i+3)+da*dx(i+3)
1140 END DO
1141END SUBROUTINE dbaxpy
1142
1151
1152SUBROUTINE dbsvx(v,a,b,n) !
1153 USE mpdef
1154
1155 IMPLICIT NONE
1156 INTEGER(mpi) :: i
1157 INTEGER(mpi) :: ij
1158 INTEGER(mpi) :: ijs
1159 INTEGER(mpi) :: j
1160
1161 ! B := V * A
1162 ! N N*N N
1163
1164 INTEGER(mpi), INTENT(IN) :: n
1165 REAL(mpd), INTENT(IN) :: v(*)
1166 REAL(mpd), INTENT(IN) :: a(*)
1167 REAL(mpd), INTENT(OUT) :: b(*)
1168
1169 REAL(mpd) :: dsum
1170 ijs=1
1171 DO i=1,n
1172 dsum=0.0
1173 ij=ijs
1174 DO j=1,n
1175 dsum=dsum+v(ij)*a(j)
1176 IF(j < i) THEN
1177 ij=ij+1
1178 ELSE
1179 ij=ij+j
1180 END IF
1181 END DO
1182 b(i)=dsum
1183 ijs=ijs+i
1184 END DO
1185END SUBROUTINE dbsvx
1186
1195
1196SUBROUTINE dbsvxl(v,a,b,n) ! LARGE symm. matrix, vector
1197 USE mpdef
1198
1199 IMPLICIT NONE
1200 INTEGER(mpi) :: i
1201 INTEGER(mpi) :: j
1202
1203 ! B := V * A
1204 ! N N*N N
1205
1206 INTEGER(mpi), INTENT(IN) :: n
1207 REAL(mpd), INTENT(IN) :: v(*)
1208 REAL(mpd), INTENT(IN) :: a(*)
1209 REAL(mpd), INTENT(OUT) :: b(*)
1210
1211 REAL(mpd) :: dsum
1212 INTEGER(mpl) :: ij
1213 INTEGER(mpl) :: ijs
1214 ijs=1
1215 DO i=1,n
1216 dsum=0.0
1217 ij=ijs
1218 DO j=1,n
1219 dsum=dsum+v(ij)*a(j)
1220 IF(j < i) THEN
1221 ij=ij+1
1222 ELSE
1223 ij=ij+int8(j)
1224 END IF
1225 END DO
1226 b(i)=dsum
1227 ijs=ijs+int8(i)
1228 END DO
1229END SUBROUTINE dbsvxl
1230
1238
1239SUBROUTINE dbgax(a,x,y,m,n)
1240 USE mpdef
1241
1242 IMPLICIT NONE
1243 INTEGER(mpi) :: i
1244 INTEGER(mpi) :: ij
1245 INTEGER(mpi) :: j
1246
1247 REAL(mpd), INTENT(IN) :: a(*)
1248 REAL(mpd), INTENT(IN) :: x(*)
1249 REAL(mpd), INTENT(OUT) :: y(*)
1250 INTEGER(mpi), INTENT(IN) :: m
1251 INTEGER(mpi), INTENT(IN) :: n
1252
1253 ! ...
1254 ij=0
1255 DO i=1,m
1256 y(i)=0.0_mpd
1257 DO j=1,n
1258 ij=ij+1
1259 y(i)=y(i)+a(ij)*x(j)
1260 END DO
1261 END DO
1262END SUBROUTINE dbgax
1263
1276SUBROUTINE dbavat(v,a,w,n,ms)
1277 USE mpdef
1278
1279 IMPLICIT NONE
1280 INTEGER(mpi) :: i
1281 INTEGER(mpi) :: ij
1282 INTEGER(mpi) :: ijs
1283 INTEGER(mpi) :: il
1284 INTEGER(mpi) :: j
1285 INTEGER(mpi) :: jk
1286 INTEGER(mpi) :: k
1287 INTEGER(mpi) :: l
1288 INTEGER(mpi) :: lk
1289 INTEGER(mpi) :: lkl
1290 INTEGER(mpi) :: m
1291
1292 REAL(mpd), INTENT(IN) :: v(*)
1293 REAL(mpd), INTENT(IN) :: a(*)
1294 REAL(mpd), INTENT(INOUT) :: w(*)
1295 INTEGER(mpi), INTENT(IN) :: n
1296 INTEGER(mpi), INTENT(IN) :: ms
1297
1298 REAL(mpd) :: cik
1299 ! ...
1300 m=ms
1301 IF (m > 0) THEN
1302 DO i=1,(m*m+m)/2
1303 w(i)=0.0_mpd ! reset output matrix
1304 END DO
1305 ELSE
1306 m=-m
1307 END IF
1308
1309 il=-n
1310 ijs=0
1311 DO i=1,m ! do I
1312 ijs=ijs+i-1 !
1313 il=il+n !
1314 lkl=0 !
1315 DO k=1,n ! do K
1316 cik=0.0_mpd !
1317 lkl=lkl+k-1 !
1318 lk=lkl !
1319 DO l=1,k ! do L
1320 lk=lk+1 ! .
1321 cik=cik+a(il+l)*v(lk) ! .
1322 END DO ! end do L
1323 DO l=k+1,n ! do L
1324 lk=lk+l-1 ! .
1325 cik=cik+a(il+l)*v(lk) ! .
1326 END DO ! end do L
1327 jk=k !
1328 ij=ijs !
1329 DO j=1,i ! do J
1330 ij=ij+1 ! .
1331 w(ij)=w(ij)+cik*a(jk) ! .
1332 jk=jk+n ! .
1333 END DO ! end do J
1334 END DO ! end do K
1335 END DO ! end do I
1336END SUBROUTINE dbavat
1337
1347
1348SUBROUTINE dbmprv(lun,x,v,n)
1349 USE mpdef
1350
1351 IMPLICIT NONE
1352 INTEGER(mpi) :: i
1353 INTEGER(mpi) :: ii
1354 INTEGER(mpi) :: ij
1355 INTEGER(mpi) :: j
1356 INTEGER(mpi) :: jj
1357 INTEGER(mpi) :: l
1358 INTEGER(mpi) :: m
1359 INTEGER(mpi) :: mc(15)
1360 REAL(mps) :: pd
1361 REAL(mps) :: rho
1362 REAL(mps) :: err
1363
1364 INTEGER(mpi), INTENT(IN) :: lun
1365 REAL(mpd), INTENT(IN) :: x(*)
1366 REAL(mpd), INTENT(IN) :: v(*)
1367 INTEGER(mpi), INTENT(IN) :: n
1368
1369 WRITE(lun,103)
1370 WRITE(lun,101)
1371 ii=0
1372 DO i=1,n
1373 ij=ii
1374 ii=ii+i
1375 err=0.0
1376 IF(v(ii) > 0.0) err=sqrt(real(v(ii),mps))
1377 l=0
1378 jj=0
1379 DO j=1,i
1380 jj=jj+j
1381 ij=ij+1
1382 rho=0.0
1383 pd=real(v(ii)*v(jj),mps)
1384 IF(pd > 0.0) rho=real(v(ij),mps)/sqrt(pd)
1385 l=l+1
1386 mc(l)=nint(100.0*abs(rho),mpi)
1387 IF(rho < 0.0) mc(l)=-mc(l)
1388 IF(j == i.OR.l == 15) THEN
1389 IF(j <= 15) THEN
1390 IF(j == i) THEN
1391 WRITE(lun,102) i,x(i),err,(mc(m),m=1,l-1)
1392 ELSE
1393 WRITE(lun,102) i,x(i),err,(mc(m),m=1,l)
1394 END IF
1395 ELSE
1396 IF(j == i) THEN
1397 WRITE(lun,103) (mc(m),m=1,l-1)
1398 ELSE
1399 WRITE(lun,103) (mc(m),m=1,l)
1400 END IF
1401 l=0
1402 END IF
1403 END IF
1404 END DO
1405 END DO
1406 WRITE(lun,104)
1407 ! 100 RETURN
1408 RETURN
1409101 FORMAT(9x,'Param',7x,'error',7x,'correlation coefficients'/)
1410102 FORMAT(1x,i5,2g12.4,1x,15i5)
1411103 FORMAT(31x,15i5)
1412104 FORMAT(33x,'(correlation coefficients in percent)')
1413END SUBROUTINE dbmprv
1414
1422
1423SUBROUTINE dbprv(lun,v,n)
1424 USE mpdef
1425
1426 IMPLICIT NONE
1427 INTEGER(mpi), PARAMETER :: istp=6
1428 INTEGER(mpi) :: i
1429 INTEGER(mpi) :: ip
1430 INTEGER(mpi) :: ipe
1431 INTEGER(mpi) :: ipn
1432 INTEGER(mpi) :: ips
1433 INTEGER(mpi) :: k
1434
1435 INTEGER(mpi), INTENT(IN) :: lun
1436 REAL(mpd), INTENT(IN) :: v(*)
1437 INTEGER(mpi), INTENT(IN) :: n
1438
1439 WRITE(lun,101)
1440
1441 DO i=1,n
1442 ips=(i*i-i)/2
1443 ipe=ips+i
1444 ip =ips
1445100 CONTINUE
1446 ipn=ip+istp
1447 WRITE(lun,102) i, ip+1-ips, (v(k),k=ip+1,min(ipn,ipe))
1448 IF (ipn < ipe) THEN
1449 ip=ipn
1450 GO TO 100
1451 END IF
1452END DO
1453RETURN
1454101 FORMAT(1x,'--- DBPRV -----------------------------------')
1455102 FORMAT(1x,2i3,6g12.4)
1456END SUBROUTINE dbprv
1457
1458! sort -------------------------------------------------------------
1459
1466
1467SUBROUTINE heapf(a,n)
1468 USE mpdef
1469
1470 IMPLICIT NONE
1471 !
1472 INTEGER(mpi) :: i
1473 INTEGER(mpi) :: j
1474 INTEGER(mpi) :: l
1475 INTEGER(mpi) :: r
1476 REAL(mps) :: at ! pivot key value
1477
1478 REAL(mps), INTENT(IN OUT) :: a(*)
1479 INTEGER(mpi), INTENT(IN) :: n
1480 ! ...
1481 IF(n <= 1) RETURN
1482 l=n/2+1
1483 r=n
148410 IF(l > 1) THEN
1485 l=l-1
1486 at =a(l)
1487 ELSE
1488 at =a(r)
1489 a(r)=a(1)
1490 r=r-1
1491 IF(r == 1) THEN
1492 a(1)=at
1493 RETURN
1494 END IF
1495 END IF
1496 i=l
1497 j=l+l
149820 IF(j <= r) THEN
1499 IF(j < r) THEN
1500 IF(a(j) < a(j+1)) j=j+1
1501 END IF
1502 IF(at < a(j)) THEN
1503 a(i)=a(j)
1504 i=j
1505 j=j+j
1506 ELSE
1507 j=r+1
1508 END IF
1509 GO TO 20
1510 END IF
1511 a(i)=at
1512 GO TO 10
1513END SUBROUTINE heapf
1514
1521
1522SUBROUTINE sort1k(a,n)
1523 USE mpdef
1524
1525 IMPLICIT NONE
1526 INTEGER(mpi) :: nlev ! stack size
1527 parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1528 INTEGER(mpi) :: i
1529 INTEGER(mpi) :: j
1530 INTEGER(mpi) :: l
1531 INTEGER(mpi) :: r
1532 INTEGER(mpi) :: lev
1533 INTEGER(mpi) :: lr(nlev)
1534 INTEGER(mpi) :: lrh
1535 INTEGER(mpi) :: maxlev
1536 INTEGER(mpi) :: a1 ! pivot key
1537 INTEGER(mpi) :: at ! pivot key
1538
1539 INTEGER(mpi), INTENT(IN OUT) :: a(*)
1540 INTEGER(mpi), INTENT(IN) :: n
1541 ! ...
1542 IF (n <= 0) RETURN
1543 maxlev=0
1544 lev=0
1545 l=1
1546 r=n
154710 IF(r-l == 1) THEN ! sort two elements L and R
1548 IF(a(l) > a(r)) THEN
1549 at=a(l) ! exchange L <-> R
1550 a(l)=a(r)
1551 a(r)=at
1552 END IF
1553 r=l
1554 END IF
1555 IF(r == l) THEN
1556 IF(lev <= 0) THEN
1557 ! WRITE(*,*) 'SORT1K (quicksort): maxlevel used/available =',
1558 ! + MAXLEV,'/64'
1559 RETURN
1560 END IF
1561 lev=lev-2
1562 l=lr(lev+1)
1563 r=lr(lev+2)
1564 ELSE
1565 ! LRH=(L+R)/2
1566 lrh=(l/2)+(r/2) ! avoid bit overflow
1567 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1568 a1=a(lrh) ! middle
1569 i=l-1 ! find limits [J,I] with [L,R]
1570 j=r+1
157120 i=i+1
1572 IF(a(i) < a1) GO TO 20
157330 j=j-1
1574 IF(a(j) > a1) GO TO 30
1575 IF(i <= j) THEN
1576 at=a(i) ! exchange I <-> J
1577 a(i)=a(j)
1578 a(j)=at
1579 GO TO 20
1580 END IF
1581 IF(lev+2 > nlev) THEN
1582 CALL peend(33,'Aborted, stack overflow in quicksort')
1583 stop 'SORT1K (quicksort): stack overflow'
1584 END IF
1585 IF(r-i < j-l) THEN
1586 lr(lev+1)=l
1587 lr(lev+2)=j
1588 l=i
1589 ELSE
1590 lr(lev+1)=i
1591 lr(lev+2)=r
1592 r=j
1593 END IF
1594 lev=lev+2
1595 maxlev=max(maxlev,lev)
1596 END IF
1597 GO TO 10
1598END SUBROUTINE sort1k
1599
1606
1607SUBROUTINE sort2k(a,n)
1608 USE mpdef
1609
1610 IMPLICIT NONE
1611 INTEGER(mpi) :: nlev ! stack size
1612 parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1613 INTEGER(mpi) :: i
1614 INTEGER(mpi) ::j
1615 INTEGER(mpi) ::l
1616 INTEGER(mpi) ::r
1617 INTEGER(mpi) ::lev
1618 INTEGER(mpi) ::lr(nlev)
1619 INTEGER(mpi) ::lrh
1620 INTEGER(mpi) ::maxlev
1621 INTEGER(mpi) ::a1 ! pivot key
1622 INTEGER(mpi) ::a2 ! pivot key
1623 INTEGER(mpi) ::at ! pivot key
1624
1625 INTEGER(mpi), INTENT(IN OUT) :: a(2,*)
1626 INTEGER(mpi), INTENT(IN) :: n
1627 ! ...
1628 maxlev=0
1629 lev=0
1630 l=1
1631 r=n
163210 IF(r-l == 1) THEN ! sort two elements L and R
1633 IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
1634 at=a(1,l) ! exchange L <-> R
1635 a(1,l)=a(1,r)
1636 a(1,r)=at
1637 at=a(2,l)
1638 a(2,l)=a(2,r)
1639 a(2,r)=at
1640 END IF
1641 r=l
1642 END IF
1643 IF(r == l) THEN
1644 IF(lev <= 0) THEN
1645 WRITE(*,*) 'SORT2K (quicksort): maxlevel used/available =', maxlev,'/64'
1646 RETURN
1647 END IF
1648 lev=lev-2
1649 l=lr(lev+1)
1650 r=lr(lev+2)
1651 ELSE
1652 ! LRH=(L+R)/2
1653 lrh=(l/2)+(r/2) ! avoid bit overflow
1654 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1655 a1=a(1,lrh) ! middle
1656 a2=a(2,lrh)
1657 i=l-1 ! find limits [J,I] with [L,R]
1658 j=r+1
165920 i=i+1
1660 IF(a(1,i) < a1) GO TO 20
1661 IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20
166230 j=j-1
1663 IF(a(1,j) > a1) GO TO 30
1664 IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
1665 IF(i <= j) THEN
1666 at=a(1,i) ! exchange I <-> J
1667 a(1,i)=a(1,j)
1668 a(1,j)=at
1669 at=a(2,i)
1670 a(2,i)=a(2,j)
1671 a(2,j)=at
1672 GO TO 20
1673 END IF
1674 IF(lev+2 > nlev) THEN
1675 CALL peend(33,'Aborted, stack overflow in quicksort')
1676 stop 'SORT2K (quicksort): stack overflow'
1677 END IF
1678 IF(r-i < j-l) THEN
1679 lr(lev+1)=l
1680 lr(lev+2)=j
1681 l=i
1682 ELSE
1683 lr(lev+1)=i
1684 lr(lev+2)=r
1685 r=j
1686 END IF
1687 lev=lev+2
1688 maxlev=max(maxlev,lev)
1689 END IF
1690 GO TO 10
1691END SUBROUTINE sort2k
1692
1699
1700SUBROUTINE sort2i(a,n)
1701 USE mpdef
1702
1703 IMPLICIT NONE
1704 INTEGER(mpi) :: nlev ! stack size
1705 parameter(nlev=2*32) ! ... for N = 2**32 = 4.3 10**9
1706 INTEGER(mpi) :: i
1707 INTEGER(mpi) ::j
1708 INTEGER(mpi) ::l
1709 INTEGER(mpi) ::r
1710 INTEGER(mpi) ::lev
1711 INTEGER(mpi) ::lr(nlev)
1712 INTEGER(mpi) ::lrh
1713 INTEGER(mpi) ::maxlev
1714 INTEGER(mpi) ::a1 ! pivot key
1715 INTEGER(mpi) ::a2 ! pivot key
1716 INTEGER(mpi) ::at ! pivot key
1717
1718 INTEGER(mpi), INTENT(IN OUT) :: a(3,*)
1719 INTEGER(mpi), INTENT(IN) :: n
1720 ! ...
1721 maxlev=0
1722 lev=0
1723 l=1
1724 r=n
172510 IF(r-l == 1) THEN ! sort two elements L and R
1726 IF(a(1,l) > a(1,r).OR.( a(1,l) == a(1,r).AND.a(2,l) > a(2,r))) THEN
1727 at=a(1,l) ! exchange L <-> R
1728 a(1,l)=a(1,r)
1729 a(1,r)=at
1730 at=a(2,l)
1731 a(2,l)=a(2,r)
1732 a(2,r)=at
1733 at=a(3,l)
1734 a(3,l)=a(3,r)
1735 a(3,r)=at
1736 END IF
1737 r=l
1738 END IF
1739 IF(r == l) THEN
1740 IF(lev <= 0) THEN
1741 WRITE(*,*) 'SORT2I (quicksort): maxlevel used/available =', maxlev,'/64'
1742 RETURN
1743 END IF
1744 lev=lev-2
1745 l=lr(lev+1)
1746 r=lr(lev+2)
1747 ELSE
1748 ! LRH=(L+R)/2
1749 lrh=(l/2)+(r/2) ! avoid bit overflow
1750 IF(mod(l,2) == 1.AND.mod(r,2) == 1) lrh=lrh+1
1751 a1=a(1,lrh) ! middle
1752 a2=a(2,lrh)
1753 i=l-1 ! find limits [J,I] with [L,R]
1754 j=r+1
175520 i=i+1
1756 IF(a(1,i) < a1) GO TO 20
1757 IF(a(1,i) == a1.AND.a(2,i) < a2) GO TO 20
175830 j=j-1
1759 IF(a(1,j) > a1) GO TO 30
1760 IF(a(1,j) == a1.AND.a(2,j) > a2) GO TO 30
1761 IF(i <= j) THEN
1762 at=a(1,i) ! exchange I <-> J
1763 a(1,i)=a(1,j)
1764 a(1,j)=at
1765 at=a(2,i)
1766 a(2,i)=a(2,j)
1767 a(2,j)=at
1768 at=a(3,i)
1769 a(3,i)=a(3,j)
1770 a(3,j)=at
1771 GO TO 20
1772 END IF
1773 IF(lev+2 > nlev) THEN
1774 CALL peend(33,'Aborted, stack overflow in quicksort')
1775 stop 'SORT2I (quicksort): stack overflow'
1776 END IF
1777 IF(r-i < j-l) THEN
1778 lr(lev+1)=l
1779 lr(lev+2)=j
1780 l=i
1781 ELSE
1782 lr(lev+1)=i
1783 lr(lev+2)=r
1784 r=j
1785 END IF
1786 lev=lev+2
1787 maxlev=max(maxlev,lev)
1788 END IF
1789 GO TO 10
1790END SUBROUTINE sort2i
1791
1799
1800REAL(mps) FUNCTION chindl(n,nd)
1801 USE mpdef
1802
1803 IMPLICIT NONE
1804 INTEGER(mpi) :: m
1805 INTEGER(mpi), INTENT(IN) :: n
1806 INTEGER(mpi), INTENT(IN) :: nd
1807 !
1808 REAL(mps) :: sn(3)
1809 REAL(mps) ::table(30,3)
1810 ! REAL PN(3)
1811 ! DATA PN/0.31731,0.0455002785,2.69985E-3/ ! probabilities
1812 DATA sn/0.47523,1.690140,2.782170/
1813 DATA table/ 1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630, &
1814 1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321, &
1815 1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136, &
1816 1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040, &
1817 4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736, &
1818 1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658, &
1819 1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346, &
1820 1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742, &
1821 9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468, &
1822 2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635, &
1823 2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901, &
1824 1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681/
1825 SAVE sn,table
1826 ! ...
1827 IF(nd < 1) THEN
1828 chindl=0.0
1829 ELSE
1830 m=max(1,min(n,3)) ! 1, 2 or 3 sigmas
1831 IF(nd <= 30) THEN
1832 chindl=table(nd,m) ! from table
1833 ELSE ! approximation for ND > 30
1834 chindl=(sn(m)+sqrt(real(nd+nd-1,mps)))**2/real(nd+nd,mps)
1835 END IF
1836 END IF
1837END FUNCTION chindl
1838
1876
1877SUBROUTINE lltdec(n,c,india,nrkd,iopt)
1878 USE mpdef
1879
1880 IMPLICIT NONE
1881 INTEGER(mpi) :: i
1882 INTEGER(mpi) :: j
1883 INTEGER(mpi) :: k
1884 INTEGER(mpi) :: kj
1885 INTEGER(mpi) :: mj
1886 INTEGER(mpi) :: mk
1887 REAL(mpd) ::diag
1888
1889 INTEGER(mpi), INTENT(IN) :: n
1890 REAL(mpd), INTENT(IN OUT) :: c(*)
1891 INTEGER(mpi), INTENT(IN) :: india(n)
1892 INTEGER(mpi), INTENT(OUT) :: nrkd
1893 INTEGER(mpi), INTENT(IN) :: iopt
1894 REAL(mpd) eps
1895 ! ...
1896 eps = 16.0_mpd * epsilon(eps) ! 16 * precision(mpd)
1897
1898 ! ..
1899 nrkd=0
1900 diag=0.0_mpd
1901 IF(c(india(1)) > 0.0) THEN
1902 c(india(1))=1.0_mpd/sqrt(c(india(1))) ! square root
1903 ELSE
1904 c(india(1))=0.0_mpd
1905 nrkd=-1
1906 END IF
1907
1908 DO k=2,n
1909 mk=k-india(k)+india(k-1)+1 ! first index in row K
1910 DO j=mk,k ! loop over row K with index J
1911 mj=1
1912 IF(j > 1) mj=j-india(j)+india(j-1)+1 ! first index in row J
1913 kj=india(k)-k+j ! index kj
1914 diag=c(india(j)) ! j-th diagonal element
1915
1916 DO i=max(mj,mk),j-1
1917 ! L_kj = L_kj - L_ki *D_ii *L_ji
1918 c(kj)=c(kj) - c(india(k)-k+i)*c(india(j)-j+i)
1919 END DO ! I
1920
1921 IF(j /= k) c(kj)=c(kj)*diag
1922 END DO ! J
1923
1924 IF(c(india(k)) > eps*diag) THEN ! test for linear dependence
1925 c(india(k))=1.0_mpd/sqrt(c(india(k))) ! square root
1926 ELSE
1927 DO j=mk,k ! reset row K
1928 c(india(k)-k+j)=0.0_mpd
1929 END DO ! J
1930 IF (iopt > 0 .and. diag > 0.0) THEN ! skyline
1931 c(india(k))=1.0_mpd/sqrt(diag) ! square root
1932 ELSE
1933 IF(nrkd == 0) THEN
1934 nrkd=-k
1935 ELSE
1936 IF(nrkd < 0) nrkd=1
1937 nrkd=nrkd+1
1938 END IF
1939 END IF
1940 END IF
1941
1942 END DO ! K
1943 RETURN
1944END SUBROUTINE lltdec
1945
1946
1958
1959SUBROUTINE lltfwd(n,c,india,x)
1960 USE mpdef
1961
1962 IMPLICIT NONE
1963 INTEGER(mpi) :: j
1964 INTEGER(mpi) :: k
1965
1966 INTEGER(mpi), INTENT(IN) :: n
1967 REAL(mpd), INTENT(IN) :: c(*)
1968 INTEGER(mpi), INTENT(IN) :: india(n)
1969 REAL(mpd), INTENT(IN OUT) :: x(n)
1970
1971 x(1)=x(1)*c(india(1))
1972 DO k=2,n ! forward loop
1973 DO j=k-india(k)+india(k-1)+1,k-1
1974 x(k)=x(k)-c(india(k)-k+j)*x(j) ! X_k := X_k - L_kj * B_j
1975 END DO ! J
1976 x(k)=x(k)*c(india(k))
1977 END DO ! K
1978
1979 RETURN
1980END SUBROUTINE lltfwd
1981
1993
1994SUBROUTINE lltbwd(n,c,india,x)
1995 USE mpdef
1996
1997 IMPLICIT NONE
1998 INTEGER(mpi) :: j
1999 INTEGER(mpi) :: k
2000
2001 INTEGER(mpi), INTENT(IN) :: n
2002 REAL(mpd), INTENT(IN) :: c(*)
2003 INTEGER(mpi), INTENT(IN) :: india(n)
2004 REAL(mpd), INTENT(IN OUT) :: x(n)
2005
2006 DO k=n,2,-1 ! backward loop
2007 x(k)=x(k)*c(india(k))
2008 DO j=k-india(k)+india(k-1)+1,k-1
2009 x(j)=x(j)-c(india(k)-k+j)*x(k) ! X_j := X_j - L_kj * X_k
2010 END DO ! J
2011 END DO ! K
2012 x(1)=x(1)*c(india(1))
2013
2014 RETURN
2015END SUBROUTINE lltbwd
2016
2034SUBROUTINE equdec(n,m,ls,c,india,nrkd,nrkd2)
2035 USE mpdef
2036
2037 IMPLICIT NONE
2038 INTEGER(mpi) :: i
2039 INTEGER(mpi) :: j
2040 INTEGER(mpi) :: jk
2041 INTEGER(mpi) :: k
2042
2043 INTEGER(mpi), INTENT(IN) :: n
2044 INTEGER(mpi), INTENT(IN) :: m
2045 INTEGER(mpi), INTENT(IN) :: ls
2046 REAL(mpd), INTENT(IN OUT) :: c(*)
2047 INTEGER(mpi), INTENT(IN OUT) :: india(n+m)
2048 INTEGER(mpi), INTENT(OUT) :: nrkd
2049 INTEGER(mpi), INTENT(OUT) :: nrkd2
2050
2051 ! ...
2052
2053 nrkd=0
2054 nrkd2=0
2055
2056 CALL lltdec(n,c,india,nrkd,ls) ! decomposition G G^T
2057
2058 IF (m>0) THEN
2059 DO i=1,m
2060 CALL lltfwd(n,c,india,c(india(n)+(i-1)*n+1)) ! forward solution K
2061 END DO
2062
2063 jk=india(n)+n*m
2064 DO j=1,m
2065 DO k=1,j
2066 jk=jk+1
2067 c(jk)=0.0_mpd ! product K K^T
2068 DO i=1,n
2069 c(jk)=c(jk)+c(india(n)+(j-1)*n+i)*c(india(n)+(k-1)*n+i)
2070 END DO
2071 END DO
2072 END DO
2073
2074 india(n+1)=1
2075 DO i=2,m
2076 india(n+i)=india(n+i-1)+min(i,m) ! pointer for K K^T
2077 END DO
2078
2079 CALL lltdec(m,c(india(n)+n*m+1),india(n+1),nrkd2,0) ! decomp. H H^T
2080 ENDIF
2081
2082 RETURN
2083END SUBROUTINE equdec
2084
2100SUBROUTINE equslv(n,m,c,india,x) ! solution vector
2101 USE mpdef
2102
2103 IMPLICIT NONE
2104 INTEGER(mpi) :: i
2105 INTEGER(mpi) :: j
2106
2107 INTEGER(mpi), INTENT(IN) :: n
2108 INTEGER(mpi), INTENT(IN) :: m
2109 REAL(mpd), INTENT(IN) :: c(*)
2110 INTEGER(mpi), INTENT(IN) :: india(n+m)
2111 REAL(mpd), INTENT(IN OUT) :: x(n+m)
2112
2113 CALL lltfwd(n,c,india,x) ! result is u
2114
2115 IF (m>0) THEN
2116 DO i=1,m
2117 DO j=1,n
2118 x(n+i)=x(n+i)-x(j)*c(india(n)+(i-1)*n+j) ! g - K u
2119 END DO
2120 END DO
2121 CALL lltfwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is v
2122
2123
2124 CALL lltbwd(m,c(india(n)+n*m+1),india(n+1),x(n+1)) ! result is -y
2125 DO i=1,m
2126 x(n+i)=-x(n+i) ! result is +y
2127 END DO
2128
2129 DO i=1,n
2130 DO j=1,m
2131 x(i)=x(i)-x(n+j)*c(india(n)+(j-1)*n+i) ! u - K^T y
2132 END DO
2133 END DO
2134 ENDIF
2135
2136 CALL lltbwd(n,c,india,x) ! result is x
2137
2138 RETURN
2139END SUBROUTINE equslv
2140
2170
2171SUBROUTINE precon(p,n,c,cu,a,s,nrkd)
2172 USE mpdef
2173
2174 IMPLICIT NONE
2175 INTEGER(mpi) :: i
2176 INTEGER(mpi) :: ii
2177 INTEGER(mpi) :: j
2178 INTEGER(mpi) :: jj
2179 INTEGER(mpi) :: jk
2180 INTEGER(mpi) :: k
2181 INTEGER(mpi) :: kk
2182
2183 INTEGER(mpi), INTENT(IN) :: p
2184 INTEGER(mpi), INTENT(IN) :: n
2185 REAL(mpd), INTENT(IN) :: c(n)
2186 REAL(mpd), INTENT(OUT) :: cu(n)
2187 REAL(mpd), INTENT(IN OUT) :: a(n,p)
2188 REAL(mpd), INTENT(OUT) :: s((p*p+p)/2)
2189 INTEGER(mpi), INTENT(OUT) :: nrkd
2190
2191 REAL(mpd) :: div
2192 REAL(mpd) :: ratio
2193
2194 nrkd=0
2195 DO i=1,(p*p+p)/2
2196 s(i)=0.0_mpd
2197 END DO
2198 DO i=1,n
2199 jk=0
2200 div=c(i) ! copy
2201 IF (div > 0.0_mpd) THEN
2202 cu(i)=1.0_mpd/sqrt(div)
2203 ELSE
2204 cu(i)=0.0_mpd
2205 nrkd=nrkd+1
2206 END IF
2207 DO j=1,p
2208 a(i,j)=a(i,j)*cu(i) ! K = A C^{-1/2}
2209 DO k=1,j
2210 jk=jk+1
2211 s(jk)=s(jk)+a(i,j)*a(i,k) ! S = symmetric matrix K K^T
2212 END DO
2213 END DO
2214 END DO
2215
2216 ii=0
2217 DO i=1,p ! S -> H D H^T (Cholesky)
2218 ii=ii+i
2219 IF(s(ii) /= 0.0_mpd) s(ii)=1.0_mpd/s(ii)
2220 jj=ii
2221 DO j=i+1,p
2222 ratio=s(i+jj)*s(ii)
2223 kk=jj
2224 DO k=j,p
2225 s(kk+j)=s(kk+j)-s(kk+i)*ratio
2226 kk=kk+k
2227 END DO ! K
2228 s(i+jj)=ratio
2229 jj=jj+j
2230 END DO ! J
2231 END DO ! I
2232 RETURN
2233END SUBROUTINE precon
2234
2244
2245SUBROUTINE presol(p,n,cu,a,s,x,y) ! solution
2246 USE mpdef
2247
2248 IMPLICIT NONE
2249 INTEGER(mpi) :: i
2250 INTEGER(mpi) :: j
2251 INTEGER(mpi) :: jj
2252 INTEGER(mpi) :: k
2253 INTEGER(mpi) :: kk
2254
2255 INTEGER(mpi), INTENT(IN) :: p
2256 INTEGER(mpi), INTENT(IN) :: n
2257
2258 REAL(mpd), INTENT(IN) :: cu(n)
2259 REAL(mpd), INTENT(IN) :: a(n,p)
2260 REAL(mpd), INTENT(IN) :: s((p*p+p)/2)
2261 REAL(mpd), INTENT(OUT) :: x(n+p)
2262 REAL(mpd), INTENT(IN) :: y(n+p)
2263
2264 REAL(mpd) :: dsum
2265
2266 DO i=1,n+p
2267 x(i)=y(i)
2268 END DO
2269 DO i=1,n
2270 x(i)=x(i)*cu(i) ! u =C^{-1/2} y
2271 DO j=1,p
2272 x(n+j)=x(n+j)-a(i,j)*x(i) ! d - K u
2273 END DO
2274 END DO
2275
2276 jj=0
2277 DO j=1,p ! Cholesky solution for v
2278 dsum=x(n+j)
2279 DO k=1,j-1
2280 dsum=dsum-s(k+jj)*x(n+k) ! H v = d - K u
2281 END DO
2282 x(n+j)=dsum ! -> v
2283 jj=jj+j
2284 END DO
2285
2286 DO j=p,1,-1 ! solution for lambda
2287 dsum=x(n+j)*s(jj)
2288 kk=jj
2289 DO k=j+1,p
2290 dsum=dsum+s(kk+j)*x(n+k) ! D H^T lambda = -v
2291 kk=kk+k
2292 END DO
2293 x(n+j)=-dsum ! -> lambda
2294 jj=jj-j
2295 END DO
2296
2297 DO i=1,n ! u - K^T lambda
2298 DO j=1,p
2299 x(i)=x(i)-a(i,j)*x(n+j)
2300 END DO
2301 END DO
2302 DO i=1,n
2303 x(i)=x(i)*cu(i) ! x = C^{-1/2} u
2304 END DO
2305
2306END SUBROUTINE presol
2307
2308
2309! 090817 C. Kleinwort, DESY-FH1
2354SUBROUTINE sqmibb(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
2355 USE mpdef
2356
2357 ! REAL(mpd) scratch arrays:
2358 ! VBND(N*(NBND+1)) = storage of band part
2359 ! VBDR(N* NBDR) = storage of border part
2360 ! AUX (N* NBDR) = intermediate results
2361
2362 ! cost[dot ops] ~= (N-NBDR)*(NBDR+NBND+1)**2 + NBDR**3/3 (leading term, solution only)
2363
2364 IMPLICIT NONE
2365 INTEGER(mpi) :: i
2366 INTEGER(mpi) :: ib
2367 INTEGER(mpi) :: ij
2368 INTEGER(mpi) :: ioff
2369 INTEGER(mpi) :: ip
2370 INTEGER(mpi) :: ip1
2371 INTEGER(mpi) :: ip2
2372 INTEGER(mpi) :: is
2373 INTEGER(mpi) :: j
2374 INTEGER(mpi) :: j0
2375 INTEGER(mpi) :: jb
2376 INTEGER(mpi) :: joff
2377 INTEGER(mpi) :: mp1
2378 INTEGER(mpi) :: nb1
2379 INTEGER(mpi) :: nmb
2380 INTEGER(mpi) :: npri
2381 INTEGER(mpi) :: nrankb
2382
2383 REAL(mpd), INTENT(IN OUT) :: v(*)
2384 REAL(mpd), INTENT(OUT) :: b(n)
2385 INTEGER(mpi), INTENT(IN) :: n
2386 INTEGER(mpi), INTENT(IN) :: nbdr
2387 INTEGER(mpi), INTENT(IN) :: nbnd
2388 INTEGER(mpi), INTENT(IN) :: inv
2389 INTEGER(mpi), INTENT(OUT) :: nrank
2390
2391 REAL(mpd), INTENT(OUT) :: vbnd(n*(nbnd+1))
2392 REAL(mpd), INTENT(OUT) :: vbdr(n*nbdr)
2393 REAL(mpd), INTENT(OUT) :: aux(n*nbdr)
2394 REAL(mpd), INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2)
2395 REAL(mpd), INTENT(OUT) :: vzru(nbdr)
2396 REAL(mpd), INTENT(OUT) :: scdiag(nbdr)
2397 INTEGER(mpi), INTENT(OUT) :: scflag(nbdr)
2398
2399 SAVE npri
2400 DATA npri / 100 /
2401 ! ...
2402 nrank=0
2403 nb1=nbdr+1
2404 mp1=nbnd+1
2405 nmb=n-nbdr
2406 ! copy band part
2407 DO i=nb1,n
2408 ip=(i*(i+1))/2
2409 is=0
2410 DO j=i,min(n,i+nbnd)
2411 ip=ip+is
2412 is=j
2413 ib=j-i+1
2414 vbnd(ib+(i-nb1)*mp1)=v(ip)
2415 END DO
2416 END DO
2417 ! copy border part
2418 IF (nbdr > 0) THEN
2419 ioff=0
2420 DO i=1,nbdr
2421 ip=(i*(i+1))/2
2422 is=0
2423 DO j=i,n
2424 ip=ip+is
2425 is=j
2426 vbdr(ioff+j)=v(ip)
2427 END DO
2428 ioff=ioff+n
2429 END DO
2430 END IF
2431
2432 CALL dbcdec(vbnd,mp1,nmb,aux)
2433 ! use? CALL DBFDEC(VBND,MP1,NMB) ! modified decomp., numerically more stable
2434 ! CALL DBCPRB(VBND,MP1,NMB)
2435 ip=1
2436 DO i=1, nmb
2437 IF (vbnd(ip) <= 0.0_mpd) THEN
2438 npri=npri-1
2439 IF (npri >= 0) THEN
2440 IF (vbnd(ip) == 0.0_mpd) THEN
2441 print *, ' SQMIBB matrix singular', n, nbdr, nbnd
2442 ELSE
2443 print *, ' SQMIBB matrix not positive definite', n, nbdr, nbnd
2444 END IF
2445 END IF
2446 ! return zeros
2447 DO ip=1,n
2448 b(ip)=0.0_mpd
2449 END DO
2450 DO ip=1,(n*n+n)/2
2451 v(ip)=0.0_mpd
2452 END DO
2453 RETURN
2454 END IF
2455 ip=ip+mp1
2456 END DO
2457 nrank=nmb
2458
2459 IF (nbdr == 0) THEN ! special case NBDR=0
2460
2461 CALL dbcslv(vbnd,mp1,nmb,b,b)
2462 IF (inv > 0) THEN
2463 IF (inv > 1) THEN
2464 CALL dbcinv(vbnd,mp1,nmb,v)
2465 ELSE
2466 CALL dbcinb(vbnd,mp1,nmb,v)
2467 END IF
2468 END IF
2469
2470 ELSE ! general case NBDR>0
2471
2472 ioff=nb1
2473 DO ib=1,nbdr
2474 ! solve for aux. vectors
2475 CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff),aux(ioff))
2476 ! zT ru
2477 vzru(ib)=b(ib)
2478 DO i=0,nmb-1
2479 vzru(ib)=vzru(ib)-b(nb1+i)*aux(ioff+i)
2480 END DO
2481 ioff=ioff+n
2482 END DO
2483 ! solve for band part only
2484 CALL dbcslv(vbnd,mp1,nmb,b(nb1),b(nb1))
2485 ! Ck - cT z
2486 ip=0
2487 ioff=nb1
2488 DO ib=1,nbdr
2489 joff=nb1
2490 DO jb=1,ib
2491 ip=ip+1
2492 vbk(ip)=v(ip)
2493 DO i=0,nmb-1
2494 vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i)
2495 END DO
2496 joff=joff+n
2497 END DO
2498 ioff=ioff+n
2499 END DO
2500 ! solve border part
2501 CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag)
2502 IF (nrankb == nbdr) THEN
2503 nrank=nrank+nbdr
2504 ELSE
2505 npri=npri-1
2506 IF (npri >= 0) print *, ' SQMIBB undef border ', n, nbdr, nbnd, nrankb
2507 DO ib=1,nbdr
2508 vzru(ib)=0.0_mpd
2509 END DO
2510 DO ip=(nbdr*nbdr+nbdr)/2,1,-1
2511 vbk(ip)=0.0_mpd
2512 END DO
2513 END IF
2514 ! smoothed data points
2515 ioff=nb1
2516 DO ib=1, nbdr
2517 b(ib) = vzru(ib)
2518 DO i=0,nmb-1
2519 b(nb1+i)=b(nb1+i)-b(ib)*aux(ioff+i)
2520 END DO
2521 ioff=ioff+n
2522 END DO
2523 ! inverse requested ?
2524 IF (inv > 0) THEN
2525 IF (inv > 1) THEN
2526 CALL dbcinv(vbnd,mp1,nmb,v)
2527 ELSE
2528 CALL dbcinb(vbnd,mp1,nmb,v)
2529 END IF
2530 ! expand/correct from NMB to N
2531 ip1=(nmb*nmb+nmb)/2
2532 ip2=(n*n+n)/2
2533 DO i=nmb-1,0,-1
2534 j0=0
2535 IF (inv == 1) j0=max(0,i-nbnd)
2536 DO j=i,j0,-1
2537 v(ip2)=v(ip1)
2538 ioff=nb1
2539 DO ib=1,nbdr
2540 joff=nb1
2541 DO jb=1,nbdr
2542 ij=max(ib,jb)
2543 ij=(ij*ij-ij)/2+min(ib,jb)
2544 v(ip2)=v(ip2)+vbk(ij)*aux(ioff+i)*aux(joff+j)
2545 joff=joff+n
2546 END DO
2547 ioff=ioff+n
2548 END DO
2549 ip1=ip1-1
2550 ip2=ip2-1
2551 END DO
2552 ip1=ip1-j0
2553 ip2=ip2-j0
2554
2555 DO ib=nbdr,1,-1
2556 v(ip2)=0.0_mpd
2557 joff=nb1
2558 DO jb=1,nbdr
2559 ij=max(ib,jb)
2560 ij=(ij*ij-ij)/2+min(ib,jb)
2561 v(ip2)=v(ip2)-vbk(ij)*aux(i+joff)
2562 joff=joff+n
2563 END DO
2564 ip2=ip2-1
2565 END DO
2566 END DO
2567
2568 DO ip=(nbdr*nbdr+nbdr)/2,1,-1
2569 v(ip2)=vbk(ip)
2570 ip2=ip2-1
2571 END DO
2572
2573 END IF
2574 END IF
2575
2576END SUBROUTINE sqmibb
2577
2578! 181105 C. Kleinwort, DESY-BELLE
2610SUBROUTINE sqmibb2(v,b,n,nbdr,nbnd,inv,nrank,vbnd,vbdr,aux,vbk,vzru,scdiag,scflag)
2611 USE mpdef
2612
2613 ! REAL(mpd) scratch arrays:
2614 ! VBND(N*(NBND+1)) = storage of band part
2615 ! VBDR(N* NBDR) = storage of border part
2616 ! AUX (N* NBDR) = intermediate results
2617
2618 ! cost[dot ops] ~= (N-NBDR)*(NBDR+NBND+1)**2 + NBDR**3/3 (leading term, solution only)
2619
2620 IMPLICIT NONE
2621 INTEGER(mpi) :: i
2622 INTEGER(mpi) :: ib
2623 INTEGER(mpi) :: ij
2624 INTEGER(mpi) :: ioff
2625 INTEGER(mpi) :: ip
2626 INTEGER(mpi) :: ip1
2627 INTEGER(mpi) :: is
2628 INTEGER(mpi) :: j
2629 INTEGER(mpi) :: j0
2630 INTEGER(mpi) :: jb
2631 INTEGER(mpi) :: joff
2632 INTEGER(mpi) :: koff
2633 INTEGER(mpi) :: mp1
2634 INTEGER(mpi) :: nb1
2635 INTEGER(mpi) :: nmb
2636 INTEGER(mpi) :: npri
2637 INTEGER(mpi) :: nrankb
2638
2639 REAL(mpd), INTENT(IN OUT) :: v(*)
2640 REAL(mpd), INTENT(OUT) :: b(n)
2641 INTEGER(mpi), INTENT(IN) :: n
2642 INTEGER(mpi), INTENT(IN) :: nbdr
2643 INTEGER(mpi), INTENT(IN) :: nbnd
2644 INTEGER(mpi), INTENT(IN) :: inv
2645 INTEGER(mpi), INTENT(OUT) :: nrank
2646
2647 REAL(mpd), INTENT(OUT) :: vbnd(n*(nbnd+1))
2648 REAL(mpd), INTENT(OUT) :: vbdr(n*nbdr)
2649 REAL(mpd), INTENT(OUT) :: aux(n*nbdr)
2650 REAL(mpd), INTENT(OUT) :: vbk((nbdr*nbdr+nbdr)/2)
2651 REAL(mpd), INTENT(OUT) :: vzru(nbdr)
2652 REAL(mpd), INTENT(OUT) :: scdiag(nbdr)
2653 INTEGER(mpi), INTENT(OUT) :: scflag(nbdr)
2654
2655 SAVE npri
2656 DATA npri / 100 /
2657 ! ...
2658 nrank=0
2659 mp1=nbnd+1
2660 nmb=n-nbdr
2661 nb1=nmb+1
2662 ! copy band part
2663 DO i=1,nmb
2664 ip=(i*(i+1))/2
2665 is=0
2666 DO j=i,min(nmb,i+nbnd)
2667 ip=ip+is
2668 is=j
2669 ib=j-i+1
2670 vbnd(ib+(i-1)*mp1)=v(ip)
2671 END DO
2672 END DO
2673 ! copy border part
2674 IF (nbdr > 0) THEN
2675 ioff=0
2676 DO i=nb1,n
2677 ip=(i*(i-1))/2
2678 DO j=1,i
2679 vbdr(ioff+j)=v(ip+j)
2680 END DO
2681 ioff=ioff+n
2682 END DO
2683 END IF
2684
2685 CALL dbcdec(vbnd,mp1,nmb,aux)
2686 ! use? CALL DBFDEC(VBND,MP1,NMB) ! modified decomp., numerically more stable
2687 ! CALL DBCPRB(VBND,MP1,NMB)
2688 ip=1
2689 DO i=1, nmb
2690 IF (vbnd(ip) <= 0.0_mpd) THEN
2691 npri=npri-1
2692 IF (npri >= 0) THEN
2693 IF (vbnd(ip) == 0.0_mpd) THEN
2694 print *, ' SQMIBB2 matrix singular', n, nbdr, nbnd
2695 ELSE
2696 print *, ' SQMIBB2 matrix not positive definite', n, nbdr, nbnd
2697 END IF
2698 END IF
2699 ! return zeros
2700 DO ip=1,n
2701 b(ip)=0.0_mpd
2702 END DO
2703 DO ip=1,(n*n+n)/2
2704 v(ip)=0.0_mpd
2705 END DO
2706 RETURN
2707 END IF
2708 ip=ip+mp1
2709 END DO
2710 nrank=nmb
2711
2712 IF (nbdr == 0) THEN ! special case NBDR=0
2713
2714 CALL dbcslv(vbnd,mp1,nmb,b,b)
2715 IF (inv > 0) THEN
2716 IF (inv > 1) THEN
2717 CALL dbcinv(vbnd,mp1,nmb,v)
2718 ELSE
2719 CALL dbcinb(vbnd,mp1,nmb,v)
2720 END IF
2721 END IF
2722
2723 ELSE ! general case NBDR>0
2724
2725 ioff=0
2726 DO ib=1,nbdr
2727 ! solve for aux. vectors
2728 CALL dbcslv(vbnd,mp1,nmb,vbdr(ioff+1),aux(ioff+1))
2729 ! zT ru
2730 vzru(ib)=b(nmb+ib)
2731 DO i=1,nmb
2732 vzru(ib)=vzru(ib)-b(i)*aux(ioff+i)
2733 END DO
2734 ioff=ioff+n
2735 END DO
2736 ! solve for band part only
2737 CALL dbcslv(vbnd,mp1,nmb,b,b)
2738 ! Ck - cT z
2739 ip=0
2740 ioff=0
2741 koff=nmb
2742 DO ib=1,nbdr
2743 joff=0
2744 DO jb=1,ib
2745 ip=ip+1
2746 vbk(ip)=vbdr(koff+jb)
2747 DO i=1,nmb
2748 vbk(ip)=vbk(ip)-vbdr(ioff+i)*aux(joff+i)
2749 END DO
2750 joff=joff+n
2751 END DO
2752 ioff=ioff+n
2753 koff=koff+n
2754 END DO
2755
2756 ! solve border part
2757 CALL sqminv(vbk,vzru,nbdr,nrankb,scdiag,scflag)
2758 IF (nrankb == nbdr) THEN
2759 nrank=nrank+nbdr
2760 ELSE
2761 npri=npri-1
2762 IF (npri >= 0) print *, ' SQMIBB2 undef border ', n, nbdr, nbnd, nrankb
2763 DO ib=1,nbdr
2764 vzru(ib)=0.0_mpd
2765 END DO
2766 DO ip=(nbdr*nbdr+nbdr)/2,1,-1
2767 vbk(ip)=0.0_mpd
2768 END DO
2769 END IF
2770 ! smoothed data points
2771 ioff=0
2772 DO ib=1, nbdr
2773 DO i=1,nmb
2774 b(i)=b(i)-vzru(ib)*aux(ioff+i)
2775 END DO
2776 ioff=ioff+n
2777 b(nmb+ib)=vzru(ib)
2778 END DO
2779 ! inverse requested ?
2780 IF (inv > 0) THEN
2781 IF (inv > 1) THEN
2782 CALL dbcinv(vbnd,mp1,nmb,v)
2783 ELSE
2784 CALL dbcinb(vbnd,mp1,nmb,v)
2785 END IF
2786 ! assemble band and border
2787 IF (nbdr > 0) THEN
2788 ! band part
2789 ip1=(nmb*nmb+nmb)/2
2790 DO i=nmb-1,0,-1
2791 j0=0
2792 IF (inv == 1) j0=max(0,i-nbnd)
2793 DO j=i,j0,-1
2794 ioff=1
2795 DO ib=1,nbdr
2796 joff=1
2797 DO jb=1,nbdr
2798 ij=max(ib,jb)
2799 ij=(ij*ij-ij)/2+min(ib,jb)
2800 v(ip1)=v(ip1)+vbk(ij)*aux(ioff+i)*aux(joff+j)
2801 joff=joff+n
2802 END DO
2803 ioff=ioff+n
2804 END DO
2805 ip1=ip1-1
2806 END DO
2807 ip1=ip1-j0
2808 END DO
2809 ! border part
2810 ip1=(nmb*nmb+nmb)/2
2811 ip=0
2812 DO ib=1,nbdr
2813 DO i=1,nmb
2814 ip1=ip1+1
2815 v(ip1)=0.0_mpd
2816 joff=0
2817 DO jb=1,nbdr
2818 ij=max(ib,jb)
2819 ij=(ij*ij-ij)/2+min(ib,jb)
2820 v(ip1)=v(ip1)-vbk(ij)*aux(i+joff)
2821 joff=joff+n
2822 END DO
2823 END DO
2824 DO jb=1,ib
2825 ip1=ip1+1
2826 ip=ip+1
2827 v(ip1)=vbk(ip)
2828 END DO
2829 END DO
2830
2831 END IF
2832 END IF
2833 END IF
2834
2835END SUBROUTINE sqmibb2
dbcdec
subroutine dbcdec(w, mp1, n, aux)
Decomposition of symmetric band matrix.
Definition Dbandmatrix.f90:420
cholin
subroutine cholin(g, v, n)
Inversion after decomposition.
Definition mpnum.f90:846
sqmibb2
subroutine sqmibb2(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Band bordered matrix.
Definition mpnum.f90:2611
vabmmm
subroutine vabmmm(n, val, ilptr)
Variable band matrix print minimum and maximum.
Definition mpnum.f90:1013
precon
subroutine precon(p, n, c, cu, a, s, nrkd)
Constrained preconditioner, decomposition.
Definition mpnum.f90:2172
devsol
subroutine devsol(n, diag, u, b, x, work)
Solution by diagonalization.
Definition mpnum.f90:658
sqminl
subroutine sqminl(v, b, n, nrank, diag, next)
Matrix inversion for LARGE matrices.
Definition mpnum.f90:225
dbsvxl
subroutine dbsvxl(v, a, b, n)
Product LARGE symmetric matrix, vector.
Definition mpnum.f90:1197
devrot
subroutine devrot(n, diag, u, v, work, iwork)
Diagonalization.
Definition mpnum.f90:378
dbmprv
subroutine dbmprv(lun, x, v, n)
Print symmetric matrix, vector.
Definition mpnum.f90:1349
dbaxpy
subroutine dbaxpy(n, da, dx, dy)
Multiply, addition.
Definition mpnum.f90:1122
equslv
subroutine equslv(n, m, c, india, x)
Solution of equilibrium systems (after decomposition).
Definition mpnum.f90:2101
heapf
subroutine heapf(a, n)
Heap sort direct (real).
Definition mpnum.f90:1468
dbdot
real(mpd) function dbdot(n, dx, dy)
Dot product.
Definition mpnum.f90:1091
equdec
subroutine equdec(n, m, ls, c, india, nrkd, nrkd2)
Decomposition of equilibrium systems.
Definition mpnum.f90:2035
vabdec
subroutine vabdec(n, val, ilptr)
Variable band matrix decomposition.
Definition mpnum.f90:901
vabslv
subroutine vabslv(n, val, ilptr, x)
Variable band matrix solution.
Definition mpnum.f90:1050
choldc
subroutine choldc(g, n)
Cholesky decomposition.
Definition mpnum.f90:754
sort1k
subroutine sort1k(a, n)
Quick sort 1.
Definition mpnum.f90:1523
sqminv
subroutine sqminv(v, b, n, nrank, diag, next)
Matrix inversion and solution.
Definition mpnum.f90:94
lltbwd
subroutine lltbwd(n, c, india, x)
Backward solution.
Definition mpnum.f90:1995
dbgax
subroutine dbgax(a, x, y, m, n)
Multiply general M-by-N matrix A and N-vector X.
Definition mpnum.f90:1240
sqmibb
subroutine sqmibb(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Bordered band matrix.
Definition mpnum.f90:2355
presol
subroutine presol(p, n, cu, a, s, x, y)
Constrained preconditioner, solution.
Definition mpnum.f90:2246
dbprv
subroutine dbprv(lun, v, n)
Print symmetric matrix.
Definition mpnum.f90:1424
devinv
subroutine devinv(n, diag, u, v)
Inversion by diagonalization. Get inverse matrix V from DIAG and U.
Definition mpnum.f90:705
lltdec
subroutine lltdec(n, c, india, nrkd, iopt)
LLT decomposition.
Definition mpnum.f90:1878
cholsl
subroutine cholsl(g, x, n)
Solution after decomposition.
Definition mpnum.f90:800
sort2k
subroutine sort2k(a, n)
Quick sort 2.
Definition mpnum.f90:1608
devsig
subroutine devsig(n, diag, u, b, coef)
Calculate significances.
Definition mpnum.f90:620
chindl
real(mps) function chindl(n, nd)
Chi2/ndf cuts.
Definition mpnum.f90:1801
lltfwd
subroutine lltfwd(n, c, india, x)
Forward solution.
Definition mpnum.f90:1960
dbsvx
subroutine dbsvx(v, a, b, n)
Product symmetric matrix, vector.
Definition mpnum.f90:1153
dbavat
subroutine dbavat(v, a, w, n, ms)
A V AT product (similarity).
Definition mpnum.f90:1277
sort2i
subroutine sort2i(a, n)
Quick sort 2 with index.
Definition mpnum.f90:1701
mpdef
Definition of constants.
Definition mpdef.f90:24
mpdef::mpd
integer, parameter mpd
Definition mpdef.f90:38
mpdef::mpi
integer, parameter mpi
Definition mpdef.f90:34
mpdef::parameter
integer, parameter parameter
Definition mpdef.f90:34
mpdef::mps
integer, parameter mps
Definition mpdef.f90:37
peend
subroutine peend(icode, cmessage)
Print exit code.
Definition pede.f90:8890