minresblas.f90

Go to the documentation of this file.

! Code converted using TO_F90 by Alan Miller
! Date: 2012-03-16  Time: 11:07:10



SUBROUTINE daxpy(n,da,dx,incx,dy,incy)

    DOUBLE PRECISION :: dx(*),dy(*),da
    INTEGER :: i,incx,incy,ix,iy,m,mp1,n

    IF(n <= 0)return
    IF (da == 0.0d0) return
    IF(incx == 1.AND.incy == 1)go to 20

    !        code for unequal increments or equal increments
    !          not equal to 1

    ix = 1
    iy = 1
    IF(incx < 0)ix = (-n+1)*incx + 1
    IF(incy < 0)iy = (-n+1)*incy + 1
    DO  i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
    END DO
    return

    !        code for both increments equal to 1


    !        clean-up loop

20  m = mod(n,4)
    IF( m == 0 ) go to 40
    DO  i = 1,m
        dy(i) = dy(i) + da*dx(i)
    END DO
    IF( n < 4 ) return
40  mp1 = m + 1
    DO  i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
    END DO

END SUBROUTINE daxpy

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


SUBROUTINE  dcopy(n,dx,incx,dy,incy)

    DOUBLE PRECISION :: dx(*),dy(*)
    INTEGER :: i,incx,incy,ix,iy,m,mp1,n

    IF(n <= 0)return
    IF(incx == 1.AND.incy == 1)go to 20

    !        code for unequal increments or equal increments
    !          not equal to 1

    ix = 1
    iy = 1
    IF(incx < 0)ix = (-n+1)*incx + 1
    IF(incy < 0)iy = (-n+1)*incy + 1
    DO  i = 1,n
        dy(iy) = dx(ix)
        ix = ix + incx
        iy = iy + incy
    END DO
    return

    !        code for both increments equal to 1


    !        clean-up loop

20  m = mod(n,7)
    IF( m == 0 ) go to 40
    DO  i = 1,m
        dy(i) = dx(i)
    END DO
    IF( n < 7 ) return
40  mp1 = m + 1
    DO  i = mp1,n,7
        dy(i) = dx(i)
        dy(i + 1) = dx(i + 1)
        dy(i + 2) = dx(i + 2)
        dy(i + 3) = dx(i + 3)
        dy(i + 4) = dx(i + 4)
        dy(i + 5) = dx(i + 5)
        dy(i + 6) = dx(i + 6)
    END DO

END SUBROUTINE dcopy

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


DOUBLE PRECISION FUNCTION ddot(n,dx,incx,dy,incy)

    DOUBLE PRECISION :: dx(*),dy(*),dtemp
    INTEGER :: i,incx,incy,ix,iy,m,mp1,n

    ddot = 0.0d0
    dtemp = 0.0d0
    IF(n <= 0)return
    IF(incx == 1.AND.incy == 1)go to 20

    !        code for unequal increments or equal increments
    !          not equal to 1

    ix = 1
    iy = 1
    IF(incx < 0)ix = (-n+1)*incx + 1
    IF(incy < 0)iy = (-n+1)*incy + 1
    DO  i = 1,n
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
    END DO
    ddot = dtemp
    return

    !        code for both increments equal to 1


    !        clean-up loop

20  m = mod(n,5)
    IF( m == 0 ) go to 40
    DO  i = 1,m
        dtemp = dtemp + dx(i)*dy(i)
    END DO
    IF( n < 5 ) go to 60
40  mp1 = m + 1
    DO  i = mp1,n,5
        dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +  &
        dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
    END DO
60  ddot = dtemp

END FUNCTION ddot

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


DOUBLE PRECISION   FUNCTION dnrm2 ( n, x, incx )

    implicit           DOUBLE PRECISION (a-h,o-z)
    INTEGER :: incx, n
    DOUBLE PRECISION :: x(*)

    !!!   double precision   s1flmx
    parameter(one = 1.0d+0, zero = 0.0d+0 )
    DOUBLE PRECISION :: norm
    INTRINSIC          abs
    !     ------------------------------------------------------------------
    !     flmax = s1flmx( )
    flmax = 1.0d+50

    IF (     n < 1) THEN
        norm  = zero
  
    ELSE IF (n == 1) THEN
        norm  = abs( x(1) )
  
    ELSE
        scale = zero
        ssq   = one
  
        DO   ix = 1, 1+(n-1)*incx, incx
    
            IF (x(ix) /= zero) THEN
                absxi = abs( x(ix) )
      
                IF (scale < absxi) THEN
                    ssq   = one + ssq*(scale/absxi)**2
                    scale = absxi
                ELSE
                    ssq   = ssq +     (absxi/scale)**2
                END IF
            END IF
        END DO

        sqt = sqrt( ssq )
        IF (scale < flmax/sqt) THEN
            norm = scale*sqt
        ELSE
            norm = flmax
        END IF
    END IF

    dnrm2  = norm

END FUNCTION dnrm2

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


SUBROUTINE daxpy2( n, a, x, y, z )

    IMPLICIT           NONE
    INTEGER :: n
    DOUBLE PRECISION :: a, x(n), y(n), z(n)
    INTEGER :: i

    DO i = 1, n
        z(i) = a*x(i) + y(i)
    END DO

END SUBROUTINE daxpy2

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


SUBROUTINE dload2( n, const, x )

    IMPLICIT           NONE
    INTEGER :: n
    DOUBLE PRECISION :: const, x(n)

    !     ------------------------------------------------------------------
    !     dload2
    !     ------------------------------------------------------------------

    INTEGER :: i

    DO i = 1, n
        x(i) = const
    END DO

END SUBROUTINE dload2

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


SUBROUTINE dscal2( n, a, x, y )

    IMPLICIT           NONE
    INTEGER :: n
    DOUBLE PRECISION :: a, x(n), y(n)

    INTEGER :: i

    DO i = 1, n
        y(i) = a*x(i)
    END DO

END SUBROUTINE dscal2