darwin3/lsopt/dgscale.F

4cee17c1be Patr*0001
                0002       subroutine dgscale( nn, gold, xdiff, diag, rmin )
                0003
                0004
                0005 c     ==================================================================
                0006 c     SUBROUTINE dgscale
                0007 c     ==================================================================
                0008 c
                0009 c     o computes new preconditioner and writes it to OPWARMD
                0010 c
                0011 c     o started: ??? not reproducible
                0012 c
                0013 c     o Version: 2.1.0, 02-Mar-2000: Patrick Heimbach, MIT/EAPS
                0014 c
                0015 c     ==================================================================
                0016 c     SUBROUTINE dgscale
                0017 c     ==================================================================
                0018
                0019       implicit none
                0020
f50e6c1777 Patr*0021 #include "blas1.h"
4cee17c1be Patr*0022
                0023       integer nn
                0024       double precision    gold(nn), xdiff(nn), diag(nn)
                0025
                0026       integer i
                0027       double precision    r1, rmin, den
                0028
                0029
                0030 c-----------------------------------------
                0031 c read diagonal
                0032 c-----------------------------------------
                0033       call dostore( nn, diag, .false., 3 )
                0034
                0035       r1 = 0.
                0036       do i = 1, nn
                0037         r1 = r1 + gold(i)*gold(i)*diag(i)
                0038       end do
                0039       r1 = 1.0 / r1
                0040
                0041       call SSCAL( nn, r1, diag, 1 )
                0042
                0043 c-----------------------------------------
                0044 c update the diagonal
                0045 c (gg is used as an auxiliary vector)
                0046 c-----------------------------------------
                0047
                0048       den = 0.0
                0049
                0050       do i = 1, nn
                0051 cph(
                0052          if (diag(i).LE.0) then
60c4531875 Patr*0053 cph            print *, 'pathei-lsopt: in dgscale; diag = 0 for i=', i
4cee17c1be Patr*0054             diag(i) = rmin
                0055          end if
                0056 cph)
                0057          den = den + xdiff(i)*xdiff(i) / diag(i)
                0058       end do
                0059
                0060       do i = 1, nn
                0061          diag(i) = 1./
                0062      $        (1./diag(i)+gold(i)**2-(xdiff(i)/diag(i))**2/den)
                0063          if (diag(i).le.0.) then
                0064             diag(i) = rmin
                0065          endif
                0066       end do
                0067
                0068 c-----------------------------------------
                0069 c write diagonal
                0070 c-----------------------------------------
                0071       call dostore( nn, diag, .true., 3 )
                0072
                0073       return
                0074       end