text:
C    Copyright (c) 2005-2017 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C        * Redistributions of source code must retain the above copyright
C          notice, this list of conditions and the following disclaimer.
C
C        * Redistributions in binary form must reproduce the above
C          copyright notice, this list of conditions and the following
C          disclaimer in the documentation and/or other materials provided
C          with the distribution.
C
C        * Neither the name of NTESS nor the names of its
C          contributors may be used to endorse or promote products derived
C          from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

      subroutine symmlq( n, b, r1, r2, v, w, x, y, work,
     $                   checka, goodb, precon, shift,
     $                   nout , itnlim, rtol, istop, itn,
     $                   anorm, acond, rnorm, ynorm, A, vwsqrt,
     $                   orthlist, macheps,normxlim,itnmin)

      external           aprod, msolve
      integer            n, nout, itnlim, istop, itn
      logical            checka, goodb, precon
      double precision   shift, rtol, anorm, acond, rnorm, ynorm,
     $                   b(n), r1(n), r2(n), v(n), w(n), x(n), y(n)
      double precision   vwsqrt(n), work(n)
      double precision   A, orthlist
      integer            itnmin
      double precision   macheps, normxlim
*     ------------------------------------------------------------------
*
*     SYMMLQ  is designed to solve the system of linear equations
*
*                Ax = b
*
*     where A is an n by n symmetric matrix and b is a given vector.
*     The matrix A is not required to be positive definite.
*     (If A is known to be definite, the method of conjugate gradients
*     might be preferred, since it will require about the same number of
*     iterations as SYMMLQ but slightly less work per iteration.)
*
*
*     The matrix A is intended to be large and sparse.  It is accessed
*     by means of a subroutine call of the form
*
*     old        call aprod ( n, x, y )
*     new:       call aprod ( n, x, y, A, vwsqrt, work, orthlist ) -rwl
*
*     which must return the product y = Ax for any given vector x.
*
*
*     More generally, SYMMLQ is designed to solve the system
*
*                (A - shift*I) x = b
*
*     where  shift  is a specified scalar value.  If  shift  and  b
*     are suitably chosen, the computed vector x may approximate an
*     (unnormalized) eigenvector of A, as in the methods of
*     inverse iteration and/or Rayleigh-quotient iteration.
*     Again, the matrix (A - shift*I) need not be positive definite.
*     The work per iteration is very slightly less if  shift = 0.
*
*
*     A further option is that of preconditioning, which may reduce
*     the number of iterations required.  If M = C C' is a positive
*     definite matrix that is known to approximate  (A - shift*I)
*     in some sense, and if systems of the form  My = x  can be
*     solved efficiently, the parameters precon and msolve may be
*     used (see below).  When  precon = .true., SYMMLQ will
*     implicitly solve the system of equations
*
*             P (A - shift*I) P' xbar  =  P b,
*
*     i.e.                  Abar xbar  =  bbar
*     where                         P  =  C**(-1),
*                                Abar  =  P (A - shift*I) P',
*                                bbar  =  P b,
*
*     and return the solution       x  =  P' xbar.
*     The associated residual is rbar  =  bbar - Abar xbar
*                                      =  P (b - (A - shift*I)x)
*                                      =  P r.
*
*     In the discussion below, eps refers to the machine precision.
*     eps is computed by SYMMLQ.  A typical value is eps = 2.22e-16
*     for IBM mainframes and IEEE double-precision arithmetic.
*
*     Parameters
*     ----------
*
*     n       input      The dimension of the matrix A.
*
*     b(n)    input      The rhs vector b.
*
*     r1(n)   workspace
*     r2(n)   workspace
*     v(n)    workspace
*     w(n)    workspace
*
*     x(n)    output     Returns the computed solution  x.
*
*     y(n)    workspace
*
*     aprod   external   A subroutine defining the matrix A.
*                        For a given vector x, the statement
*
*                        old: call aprod ( n, x, y, )
*                        new: call aprod ( n, x, y, A, vwsqrt, work, orthlist ) -rwl
*
*                        must return the product y = Ax
*                        without altering the vector x.
*
*     msolve  external   An optional subroutine defining a
*                        preconditioning matrix M, which should
*                        approximate (A - shift*I) in some sense.
*                        M must be positive definite.
*                        For a given vector x, the statement
*
*                        old:      call msolve( n, x, y )
*                        new:      call msolve( n, x, y )  -rwl
*
*                        must solve the linear system My = x
*                        without altering the vector x.
*
*                        In general, M should be chosen so that Abar has
*                        clustered eigenvalues.  For example,
*                        if A is positive definite, Abar would ideally
*                        be close to a multiple of I.
*                        If A or A - shift*I is indefinite, Abar might
*                        be close to a multiple of diag( I  -I ).
*
*                        NOTE.  The program calling SYMMLQ must declare
*                        aprod and msolve to be external.
*
*     checka  input      If checka = .true., an extra call of aprod will
*                        be used to check if A is symmetric.  Also,
*                        if precon = .true., an extra call of msolve
*                        will be used to check if M is symmetric.
*
*     goodb   input      Usually, goodb should be .false.
*                        If x is expected to contain a large multiple of
*                        b (as in Rayleigh-quotient iteration),
*                        better precision may result if goodb = .true.
*                        See Lewis (1977) below.
*                        When goodb = .true., an extra call to msolve
*                        is required.
*
*     precon  input      If precon = .true., preconditioning will
*                        be invoked.  Otherwise, subroutine msolve
*                        will not be referenced; in this case the
*                        actual parameter corresponding to msolve may
*                        be the same as that corresponding to aprod.
*
*     shift   input      Should be zero if the system Ax = b is to be
*                        solved.  Otherwise, it could be an
*                        approximation to an eigenvalue of A, such as
*                        the Rayleigh quotient b'Ab / (b'b)
*                        corresponding to the vector b.
*                        If b is sufficiently like an eigenvector
*                        corresponding to an eigenvalue near shift,
*                        then the computed x may have very large
*                        components.  When normalized, x may be
*                        closer to an eigenvector than b.
*
*     nout    input      A file number.
*                        If nout .gt. 0, a summary of the iterations
*                        will be printed on unit nout.
*
*     itnlim  input      An upper limit on the number of iterations.
*
*     rtol    input      A user-specified tolerance.  SYMMLQ terminates
*                        if it appears that norm(rbar) is smaller than
*                              rtol * norm(Abar) * norm(xbar),
*                        where rbar is the transformed residual vector,
*                              rbar = bbar - Abar xbar.
*
*                        If shift = 0 and precon = .false., SYMMLQ
*                        terminates if norm(b - A*x) is smaller than
*                              rtol * norm(A) * norm(x).
*
*     istop   output     An integer giving the reason for termination...
*
*              -1        beta2 = 0 in the Lanczos iteration; i.e. the
*                        second Lanczos vector is zero.  This means the
*                        rhs is very special.
*                        If there is no preconditioner, b is an
*                        eigenvector of A.
*                        Otherwise (if precon is true), let My = b.
*                        If shift is zero, y is a solution of the
*                        generalized eigenvalue problem Ay = lambda My,
*                        with lambda = alpha1 from the Lanczos vectors.
*
*                        In general, (A - shift*I)x = b
*                        has the solution         x = (1/alpha1) y
*                        where My = b.
*
*               0        b = 0, so the exact solution is x = 0.
*                        No iterations were performed.
*
*               1        Norm(rbar) appears to be less than
*                        the value  rtol * norm(Abar) * norm(xbar).
*                        The solution in  x  should be acceptable.
*
*               2        Norm(rbar) appears to be less than
*                        the value  eps * norm(Abar) * norm(xbar).
*                        This means that the residual is as small as
*                        seems reasonable on this machine.
*
*               3        Norm(Abar) * norm(xbar) exceeds norm(b)/eps,
*                        which should indicate that x has essentially
*                        converged to an eigenvector of A
*                        corresponding to the eigenvalue shift.
*
*               4        acond (see below) has exceeded 0.1/eps, so
*                        the matrix Abar must be very ill-conditioned.
*                        x may not contain an acceptable solution.
*
*               5        The iteration limit was reached before any of
*                        the previous criteria were satisfied.
*
*               6        The matrix defined by aprod does not appear
*                        to be symmetric.
*                        For certain vectors y = Av and r = Ay, the
*                        products y'y and r'v differ significantly.
*
*               7        The matrix defined by msolve does not appear
*                        to be symmetric.
*                        For vectors satisfying My = v and Mr = y, the
*                        products y'y and r'v differ significantly.
*
*               8        An inner product of the form  x' M**(-1) x
*                        was not positive, so the preconditioning matrix
*                        M does not appear to be positive definite.
*
*                        If istop .ge. 5, the final x may not be an
*                        acceptable solution.
*
*               9        The norm of the iterate is > than normxlim.
*                        Termination enforced on presumption that inverse
*                        iteration is being performed -rwl.
*
*     itn     output     The number of iterations performed.
*
*     anorm   output     An estimate of the norm of the matrix operator
*                        Abar = P (A - shift*I) P',   where P = C**(-1).
*
*     acond   output     An estimate of the condition of Abar above.
*                        This will usually be a substantial
*                        under-estimate of the true condition.
*
*     rnorm   output     An estimate of the norm of the final
*                        transformed residual vector,
*                           P (b  -  (A - shift*I) x).
*
*     ynorm   output     An estimate of the norm of xbar.
*                        This is sqrt( x'Mx ).  If precon is false,
*                        ynorm is an estimate of norm(x).
*
*     A       input      A pointer variable to the matrix data. Passed
*                        in to use in revised call to aprod and msolve
*                        added 14 Dec 92 by rwl.
*
*    macheps  output     Used to return the calculated machine precision.
*                        Added 10 Feb 93 by rwl.
*
*    normxlim input     Used as possible termination criterion. 10 Feb 93 rwl.
*
*    itnmin input       Used to enforce a minimum number of itns. 10 Feb 93 rwl.
*
*     To change precision
*     -------------------
*
*     Alter the words
*            double precision,
*            chdaxpy, chdcopy, chddot, chdnrm2
*     to their single or double equivalents.
*     ------------------------------------------------------------------
*
*
*     This routine is an implementation of the algorithm described in
*     the following references:
*
*     C.C. Paige and M.A. Saunders,  Solution of Sparse Indefinite
*          Systems of Linear Equations,
*          SIAM J. Numer. Anal. 12, 4, September 1975, pp. 617-629.
*
*     J.G. Lewis,  Algorithms for Sparse Matrix Eigenvalue Problems,
*          Report STAN-CS-77-595, Computer Science Department,
*          Stanford University, Stanford, California, March 1977.
*
*     Applications of SYMMLQ and the theory of preconditioning
*     are described in the following references:
*
*     D.B. Szyld and O.B. Widlund,  Applications of Conjugate Gradient
*          Type Methods to Eigenvalue Calculations,
*          in R. Vichnevetsky and R.S. Steplman (editors),
*          Advances in Computer Methods for Partial Differential
*          Equations -- III, IMACS, 1979, 167-173.
*
*     D.B. Szyld,  A Two-level Iterative Method for Large Sparse
*          Generalized Eigenvalue Calculations,
*          Ph. D. dissertation, Department of Mathematics,
*          New York University, New York, October 1983.
*
*     P.E. Gill, W. Murray, D.B. Ponceleon and M.A. Saunders,
*          Preconditioners for indefinite systems arising in
*          optimization, SIMAX 13, 1, 292--311, January 1992.
*          (SIAM J. on Matrix Analysis and Applications)
*     ------------------------------------------------------------------
*
*
*     SYMMLQ development:
*            1972: First version.
*            1975: John Lewis recommended modifications to help with
*                  inverse iteration:
*                  1. Reorthogonalize v1 and v2.
*                  2. Regard the solution as x = x1  +  bstep * b,
*                     with x1 and bstep accumulated separately
*                     and bstep * b added at the end.
*                     (In inverse iteration, b might be close to the
*                     required x already, so x1 may be a lot smaller
*                     than the multiple of b.)
*            1978: Daniel Szyld and Olof Widlund implemented the first
*                  form of preconditioning.
*                  This required both a solve and a multiply with M.
*            1979: Implemented present method for preconditioning.
*                  This requires only a solve with M.
*            1984: Sven Hammarling noted corrections to tnorm and x1lq.
*                  SYMMLQ added to NAG Fortran Library.
*     15 Sep 1985: Final F66 version.  SYMMLQ sent to "misc" in netlib.
*     16 Feb 1989: First F77 version.
*
*     22 Feb 1989: Hans Mittelmann observed beta2 = 0 (hence failure)
*                  if Abar = const*I.  istop = -1 added for this case.
*
*     01 Mar 1989: Hans Mittelmann observed premature termination on
*                  ( 1  1  1 )     (   )                   ( 1  1    )
*                  ( 1  1    ) x = ( 1 ),  for which  T3 = ( 1  1  1 ).
*                  ( 1     1 )     (   )                   (    1  1 )
*                  T2 is exactly singular, so estimating cond(A) from
*                  the diagonals of Lbar is unsafe.  We now use
*                  L       or  Lbar         depending on whether
*                  lqnorm  or  cgnorm       is least.
*
*     03 Mar 1989: eps computed internally instead of coming in as a
*                  parameter.
*     07 Jun 1989: ncheck added as a parameter to say if A and M
*                  should be checked for symmetry.
*                  Later changed to checka (see below).
*     20 Nov 1990: goodb added as a parameter to make Lewis's changes
*                  an option.  Usually b is NOT much like x.  Setting
*                  goodb = .false. saves a call to msolve at the end.
*     20 Nov 1990: Residual not computed exactly at end, to save time
*                  when only one or two iterations are required
*                  (e.g. if the preconditioner is very good).
*                  Beware, if precon is true, rnorm estimates the
*                  residual of the preconditioned system, not Ax = b.
*     04 Sep 1991: Parameter list changed and reordered.
*                  integer ncheck is now logical checka.
*     22 Jul 1992: Example from Lothar Reichel and Daniela Calvetti
*                  showed that beta2 = 0 (istop = -1) means that
*                  b is an eigenvector when M = I.
*                  More complicated if there is a preconditioner;
*                  not clear yet how to describe it.
*    14 Dec 1992:  Modified by Robert Leland, Sandia National Laboratories
*                  to integrate with a C application code. The matrix
*                  data is now passed by reference through symmlq to
*                  aprod and msolve. These are now just Fortran wrappers
*                  for C codes consistent with the matrix data passed
*                  via the pointers "A", "vwsqrt", "work" and "orthlist"
*    10 Feb 1993:  Modified by Robert Leland to return calculate machine
*                  precision and terminate if the norm of the iterate gets
*                  above the limit normxlim. Relevant for inverse iteration.
*                  Also incorporated itnmin to enforce minimum number itns.
*    17 Aug 1993:  Observed that the Fortran i/o in this routine won't
*                  work because there is no main fortran program to open
*                  the standard i/o files. So for this (and other reasons)
*                  converted the Fortran to C, necessitating inclusion of
*                  the file f2c.h. To avoid a problem with maintaining
*                  Symmlq, I commented out the i/o within it and instead
*                  report its performance based only on the return value
*                  of various parameters. That means we can modifiy the
*                  Fortran source, run f2c and recompile without losing or
*                  re-writing any functionality.
*
*     Michael A. Saunders                    na.saunders@na-net.ornl.gov
*     Department of Operations Research    mike@sol-michael.stanford.edu
*     Stanford University
*     Stanford, CA 94305-4022                             (415) 723-1875
*     ------------------------------------------------------------------
*
*
*     Subroutines and functions
*
*     USER       aprod, msolve
*     BLAS       chdaxpy, chdcopy, chddot , chdnrm2
*
*
*     Intrinsics and local variables

      intrinsic          abs, max, min, mod, sqrt
      double precision   chddot, chdnrm2
      double precision   alfa, b1, beta, beta1, bstep, cs,
     $                   cgnorm, dbar, delta, denom, diag,
     $                   eps, epsa, epsln, epsr, epsx,
     $                   gamma, gbar, gmax, gmin,
     $                   lqnorm, oldb, qrnorm, rhs1, rhs2,
     $                   s, sn, snprod, t, tnorm,
     $                   x1cg, x1lq, ynorm2, zbar, z
      integer            i

      double precision   zero         ,  one         ,  two
      parameter        ( zero = 0.0d+0,  one = 1.0d+0,  two = 2.0d+0 )

      character*16       enter, exit
      character*52       msg(-1:9)

      data               enter /' Enter SYMMLQ.  '/,
     $                   exit  /' Exit  SYMMLQ.  '/

      data               msg
     $ / 'beta2 = 0.  If M = I, b and x are eigenvectors of A',
     $   'beta1 = 0.  The exact solution is  x = 0',
     $   'Requested accuracy achieved, as determined by rtol',
     $   'Reasonable accuracy achieved, given eps',
     $   'x has converged to an eigenvector',
     $   'acond has exceeded 0.1/eps',
     $   'The iteration limit was reached',
     $   'aprod  does not define a symmetric matrix',
     $   'msolve does not define a symmetric matrix',
     $   'msolve does not define a pos-def preconditioner',
     $   'Norm of iterate > max for well conditioned system' /
*     ------------------------------------------------------------------

*     Compute eps, the machine precision.  The call to chdaxpy is
*     intended to fool compilers that use extra-length registers.

      eps    = one / 16.0d+0

   10 eps    = eps / two
      x(1)   = eps
      y(1)   = one
      call chdaxpy ( 1, one, x, 1, y, 1 )
      if (y(1) .gt. one) go to 10

      eps    = eps * two

*     Return the calculated machine precision - rwl
      macheps = eps

*     Print heading and initialize.

c   This i/o won't work - see note in preamble..
c     if (nout .gt. 0) then
c        write(nout, 1000) enter, n, checka, goodb, precon,
c     $                     itnlim, rtol, shift
c     end if

      istop  = 0
      itn    = 0
      anorm  = zero
      acond  = zero
      rnorm  = zero
      ynorm  = zero

      do 50 i = 1, n
         x(i) = zero
   50 continue

*     Set up y for the first Lanczos vector v1.
*     y is really beta1 * P * v1  where  P = C**(-1).
*     y and beta1 will be zero if b = 0.

      call chdcopy ( n, b, 1, y , 1 )
      call chdcopy ( n, b, 1, r1, 1 )
      if ( precon ) call msolve( n, r1, y, A, vwsqrt, work )
c     if ( goodb  ) then
c        b1  = y(1)
c     else
c        b1  = zero
c     end if
      beta1  = chddot  ( n, r1, 1, y, 1 )

*     See if msolve is symmetric.

      if (checka  .and.  precon) then
         call msolve( n, y, r2, A, vwsqrt, work )
         s      = chddot  ( n, y, 1, y, 1 )
         t      = chddot  ( n,r1, 1,r2, 1 )
         z      = abs( s - t )
         epsa   = (s + eps) * eps**0.33333
         if (z .gt. epsa) then
            istop = 7
            go to 900
         end if
      end if

*     Test for an indefinite preconditioner.

      if (beta1 .lt. zero) then
         istop = 8
         go to 900
      end if

*     If b = 0 exactly, stop with x = 0.

      if (beta1 .eq. zero) then
         go to 900
      end if

*     Here and later, v is really P * (the Lanczos v).

      beta1  = sqrt( beta1 )
      s      = one / beta1
      do 100 i = 1, n
         v(i)  = s * y(i)
  100 continue

*     See if aprod  is symmetric.

      call aprod( n, v, y, A, vwsqrt, work, orthlist )
      if (checka) then
         call aprod ( n, y, r2, A, vwsqrt, work, orthlist )
         s      = chddot  ( n, y, 1, y, 1 )
         t      = chddot  ( n, v, 1,r2, 1 )
         z      = abs( s - t )
         epsa   = (s + eps) * eps**0.33333
         if (z .gt. epsa) then
            istop = 6
            go to 900
         end if
      end if

*     Set up y for the second Lanczos vector.
*     Again, y is beta * P * v2  where  P = C**(-1).
*     y and beta will be zero or very small if b is an eigenvector.

      call chdaxpy ( n, (- shift), v, 1, y, 1 )
      alfa   = chddot  ( n, v, 1, y, 1 )
      call chdaxpy ( n, (- alfa / beta1), r1, 1, y, 1 )

*     Make sure  r2  will be orthogonal to the first  v.

      z      = chddot  ( n, v, 1, y, 1 )
      s      = chddot  ( n, v, 1, v, 1 )
      call chdaxpy ( n, (- z / s), v, 1, y, 1 )

      call chdcopy ( n, y, 1, r2, 1 )
      if ( precon ) call msolve( n, r2, y, A, vwsqrt, work )
      oldb   = beta1
      beta   = chddot  ( n, r2, 1, y, 1 )

      if (beta .lt. zero) then
         istop = 8
         go to 900
      end if

*     Cause termination (later) if beta is essentially zero.

      beta   = sqrt( beta )
      if (beta .le. eps) then
         istop = -1
      end if

*     See if the local reorthogonalization achieved anything.

      denom  = sqrt( s ) * chdnrm2( n, r2, 1 )  +  eps
      s      = z / denom
      t      = chddot  ( n, v, 1, r2, 1 ) / denom

c     if (nout .gt. 0  .and.  goodb) then
c        write(nout, 1100) beta1, alfa, s, t
c     end if

*     Initialize other quantities.

      cgnorm = beta1
      gbar   = alfa
      dbar   = beta
      rhs1   = beta1
      rhs2   = zero
      bstep  = zero
      snprod = one
      tnorm  = alfa**2
      ynorm2 = zero
      gmax   = abs( alfa )
      gmin   = gmax

      if ( goodb ) then
         do 200 i = 1, n
            w(i)  = zero
  200    continue
      else
         call chdcopy ( n, v, 1, w, 1 )
      end if

*     ------------------------------------------------------------------
*     Main iteration loop.
*     ------------------------------------------------------------------

*     Estimate various norms and test for convergence.

  300 anorm  = sqrt( tnorm  )
      ynorm  = sqrt( ynorm2 )
      epsa   = anorm * eps
      epsx   = anorm * ynorm * eps
      epsr   = anorm * ynorm * rtol
      diag   = gbar
      if (diag .eq. zero) diag = epsa

      lqnorm = sqrt( rhs1**2 + rhs2**2 )
      qrnorm = snprod * beta1
      cgnorm = qrnorm * beta / abs( diag )

*     Estimate  cond(A).
*     In this version we look at the diagonals of  L  in the
*     factorization of the tridiagonal matrix,  T = L*Q.
*     Sometimes, T(k) can be misleadingly ill-conditioned when
*     T(k+1) is not, so we must be careful not to overestimate acond.

      if (lqnorm .le. cgnorm) then
         acond  = gmax / gmin
      else
         denom  = min( gmin, abs( diag ) )
         acond  = gmax / denom
      end if

*     See if any of the stopping criteria are satisfied.
*     In rare cases, istop is already -1 from above (Abar = const * I).

      if (istop .eq. 0) then
         if (itn    .ge. itnlim ) istop = 5
c         if (acond  .ge. 0.1/eps) istop = 4
         if (epsx   .ge. beta1  ) istop = 3
         if (cgnorm .le. epsx   ) istop = 2
         if (cgnorm .le. epsr   ) istop = 1
         if ((istop .ne. 4).and.(ynorm .ge. normxlim)
     1        .and. (normxlim .ne. 0.0)) istop = 9
      end if
      if (itn .lt. itnmin) istop = 0
*     ==================================================================

*     See if it is time to print something.

      if (nout .le.  0)          go to 600
      if (n    .le. 40)          go to 400
      if (itn  .le. 10)          go to 400
      if (itn  .ge. itnlim - 10) go to 400
      if (mod(itn,10)  .eq.   0) go to 400
      if (cgnorm .le. 10.0*epsx) go to 400
      if (cgnorm .le. 10.0*epsr) go to 400
      if (acond  .ge. 0.01/eps ) go to 400
      if (istop  .ne. 0)         go to 400
      go to 600

*     Print a line for this iteration.

  400 zbar   = rhs1 / diag
      z      = (snprod * zbar  +  bstep) / beta1
c     x1lq   = x(1)  +  b1 * bstep / beta1
c     x1cg   = x(1)  +  w(1) * zbar  +  b1 * z

c     if (    itn     .eq. 0) write(nout, 1200)
c     write(nout, 1300) itn, x1cg, cgnorm, bstep/beta1, anorm, acond
c     if (mod(itn,10) .eq. 0) write(nout, 1500)
*     ==================================================================

*     Obtain the current Lanczos vector  v = (1 / beta)*y
*     and set up  y  for the next iteration.

  600 if (istop .ne. 0) go to 800
      s      = one / beta

      do 620 i = 1, n
         v(i)  = s * y(i)
  620 continue

      call aprod ( n, v, y, A, vwsqrt, work, orthlist )
      call chdaxpy ( n, (- shift), v, 1, y, 1 )
      call chdaxpy ( n, (- beta / oldb), r1, 1, y, 1 )
      alfa   = chddot( n, v, 1, y, 1 )
      tnorm  = tnorm  +  alfa**2  +  two * beta**2
      call chdaxpy ( n, (- alfa / beta), r2, 1, y, 1 )
      call chdcopy ( n, r2, 1, r1, 1 )
      call chdcopy ( n, y, 1, r2, 1 )
      if ( precon ) call msolve( n, r2, y, A, vwsqrt, work )
      oldb   = beta
      beta   = chddot  ( n, r2, 1, y, 1 )
      if (beta .lt. zero) then
         istop = 6
         go to 800
      end if
      beta   = sqrt( beta )

*     Compute the next plane rotation for  Q.

      gamma  = sqrt( gbar**2 + oldb**2 )
      cs     = gbar / gamma
      sn     = oldb / gamma
      delta  = cs * dbar  +  sn * alfa
      gbar   = sn * dbar  -  cs * alfa
      epsln  = sn * beta
      dbar   =            -  cs * beta

*     Update  x.

      z      = rhs1 / gamma
      s      = z * cs
      t      = z * sn

      do 700 i = 1, n
         x(i)  = (w(i) * s   +   v(i) * t)  +  x(i)
         w(i)  =  w(i) * sn  -   v(i) * cs
  700 continue

*     Accumulate the step along the direction  b,
*     and go round again.

      bstep  = snprod * cs * z  +  bstep
      snprod = snprod * sn
      gmax   = max( gmax, gamma )
      gmin   = min( gmin, gamma )
      ynorm2 = z**2  +  ynorm2
      rhs1   = rhs2  -  delta * z
      rhs2   =       -  epsln * z
      itn    = itn   +  1
      go to 300

*     ------------------------------------------------------------------
*     End of main iteration loop.
*     ------------------------------------------------------------------

*     Move to the CG point if it seems better.
*     In this version of SYMMLQ, the convergence tests involve
*     only cgnorm, so we're unlikely to stop at an LQ point,
*     EXCEPT if the iteration limit interferes.

  800 if (cgnorm .le. lqnorm) then
         zbar   = rhs1 / diag
         bstep  = snprod * zbar  +  bstep
         ynorm  = sqrt( ynorm2  +  zbar**2 )
         rnorm  = cgnorm
         call chdaxpy ( n, zbar, w, 1, x, 1 )
      else
         rnorm  = lqnorm
      end if

      if ( goodb ) then

*        Add the step along  b.

         bstep  = bstep / beta1
         call chdcopy ( n, b, 1, y, 1 )
         if ( precon ) call msolve( n, b, y, A, vwsqrt, work )
         call chdaxpy ( n, bstep, y, 1, x, 1 )
      end if

*     ==================================================================
*     Display final status.
*     ==================================================================
  900 continue
c     if (nout  .gt. 0) then
c        write(nout, 2000) exit, istop, itn,
c    $                     exit, anorm, acond,
c    $                     exit, rnorm, ynorm
c        write(nout, 3000) exit, msg(istop)
c     end if

      return

*     ------------------------------------------------------------------
c 1000 format(// 1p,    a, 5x, 'Solution of symmetric   Ax = b'
c     $       / ' n      =', i7, 5x, 'checka =', l4, 12x,
c     $          'goodb  =', l4, 7x, 'precon =', l4
c     $       / ' itnlim =', i7, 5x, 'rtol   =', e11.2, 5x,
c     $          'shift  =', e23.14)
c 1100 format(/ 1p, ' beta1  =', e10.2, 3x, 'alpha1 =', e10.2
c     $       / ' (v1,v2) before and after ', e14.2
c     $       / ' local reorthogonalization', e14.2)
c 1200 format(// 5x, 'itn', 7x, 'x1(cg)', 10x,
c     $         'norm(r)', 5x, 'bstep', 7x, 'norm(A)', 3X, 'cond(A)')
c 1300 format(1p, i8, e19.10, e11.2, e14.5, 2e10.2)
c 1500 format(1x)
c 2000 format(/ 1p, a, 6x, 'istop =', i3,   15x, 'itn   =', i8
c     $       /     a, 6x, 'anorm =', e12.4, 6x, 'acond =', e12.4
c     $       /     a, 6x, 'rnorm =', e12.4, 6x, 'ynorm =', e12.4)
c 3000 format(      a, 6x, a )
*     ------------------------------------------------------------------
*     end of SYMMLQ
      end

max_stars_repo_name: jbakosi/seacas

max_issues_repo_name: jbakosi/seacas

max_forks_repo_name: jbakosi/seacas


text:
      INTEGER*2 IFHT,IFSTRT,IFWID,IFX0,IFY0,LENGF
      BYTE NFONTS
      INTEGER IUNITF
      COMMON /PINOUT/LUNIN,LUNOUT
      COMMON /PLT$FNT/IFSTRT(150,4),LENGF(150,4),IFX0(150,4),
     +       IFY0(150,4),IFWID(150,4),IFHT(150,4),NFONTS(4,3000,4)
      DATA IUNITF/89/
C
C
      IFAIL = 0
      ITEROP = -IUNITF
      OPEN(UNIT=10,FILE='font84.ascii',STATUS='NEW')
      CALL CCPDPN (ITEROP,'font84.dat','OLD','U',80,IFAIL)
      IF (IFAIL.NE.0) THEN
        WRITE (6,FMT=6000)
 6000   FORMAT (' Unable to read fonts - FILE=PUBLIC_FONT84')
      ELSE
        READ (IUNITF) IFSTRT,LENGF,IFX0,IFY0,IFWID,IFHT,NFONTS
        WRITE(10,2000) IFSTRT,LENGF,IFX0,IFY0,IFWID,IFHT,NFONTS
2000    FORMAT(10I10)
        CLOSE (UNIT=IUNITF)
      ENDIF
      END

ext: for

max_stars_repo_name: merkys/MMB

max_issues_repo_name: merkys/MMB

max_forks_repo_name: merkys/MMB


text:
C
C
C  Copyright (C) 2000, 2001 Silicon Graphics, Inc.  All Rights Reserved.
C
C  This program is free software; you can redistribute it and/or modify it
C  under the terms of version 2.1 of the GNU Lesser General Public License
C  as published by the Free Software Foundation.
C
C  This program is distributed in the hope that it would be useful, but
C  WITHOUT ANY WARRANTY; without even the implied warranty of
C  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
C
C  Further, this software is distributed without any warranty that it is
C  free of the rightful claim of any third person regarding infringement
C  or the like.  Any license provided herein, whether implied or
C  otherwise, applies only to this software file.  Patent licenses, if
C  any, provided herein do not apply to combinations of this program with
C  other software, or any other product whatsoever.
C
C  You should have received a copy of the GNU Lesser General Public
C  License along with this program; if not, write the Free Software
C  Foundation, Inc., 59 Temple Place - Suite 330, Boston MA 02111-1307,
C  USA.
C
C  Contact information:  Silicon Graphics, Inc., 1600 Amphitheatre Pky,
C  Mountain View, CA 94043, or:
C
C  http://www.sgi.com
C
C  For further information regarding this notice, see:
C
C  http://oss.sgi.com/projects/GenInfo/NoticeExplan
C
C

       SUBROUTINE WFIND(RECALL)
CDIR$ ID "@(#) libu/waio/wfind.f	92.0	10/08/98 14:57:41"

CC     WFIND.  -  FIND  ANY BUFFERS THAT MATCH TARGET WORD ADDRESSES.
C
C      WFIND IS CALLED IF THE WORD ADDRESS DESIRED IS LESS THAN EOI
C      WORD ADDRESS.  WFINDS JOB IS TO FIND IF THE PAGES NEEDED ARE
C      STILL IN CORE.
C
C      ENTRY  -  NONE
C
C      EXIT   -  RECALL = .TRUE. IF THE CONTIGUOUS BUFFER RESIDENT
C                         INITIAL PORTION OF THE REQUEST IS AT THE END OF
C                         THE BUFFER.  THE CALLING ROUTINE CAN THEN
C                         MOVE THESE PAGES AND THEN CALL *WFIND* AGAIN
C                         TO LOOK FOR MORE PAGES ELSEWHERE IN THE
C                         BUFFER.
C
C                RECALL = .FALSE. IF ABOVE STATEMENT NOT TRUE.
C
C
C      CALLS  -  *WGETB*  TO GET BUFFERS IF NO MATCHS FOUND.
C
C                *WFLUSH* TO FLUSH BUFFERS IF ROOM AT TOP.
C
C                *READWA* TO READ FROM DISK TO BUFFERS.
C
C      METHOD:
C
C        FIRST A GLOBAL SWEEP OF ALL RESIDENT WORD ADDRESSES IS MADE.
C      IF ALL THE NECESSARY WORD ADDRESSES ARE FOUND THE SEARCH IS
C      DONE AND WE EXIT.  IF WE DONT FIND ALL THE WORD ADDRESSES WE
C      THEN SEARCH FOR CONTIGOUS TARGET WORD ADDRESSES STARTING WITH
C      THE FIRST ONE.  IF WE FIND NONE *WGETB* IS CALLED TO GET THE
C      BUFFERS AND WE THEN RETURN.  IF THE BUFFER OR BUFFERS FOUND
C      ARE AT THE END OF THE BUFFER POOL WE SET RECALL = .TRUE. AND
C      THEN WE RETURN.  IF THEY ARE NOT AT END THEN THERE IS ROOM AT
C      THE END FOR SOME OF THE PAGES THAT WE NEED TO FIND.  *WFLUSH*
C      IS CALLED TO FLUSH OUT THE BUFFERS AT THE END OF THE BUFFER
C      POOL.  IF WE ARE WRITING TO THE DISK WE NEED ONLY ENTER THE
C      WORD ADDRESSES FOR THESE NEW BUFFERS IN THE *WABP* ARRAY UNLESS
C      WE ARE DOING THE LAST BUFFER OF THIS  TRANSFER AND IT  IS A
C      PARTIAL BUFFER.  IF IT IS THEN WE MUST PULL THIS BUFFER FROM
C      DISK (ALSO THE LAST BUFFER MUST BE LESS THAN EOI) SO THAT THE
C      NEW DATA CAN BE MERGED WITH THE OLD.  IF WE ARE READING WE READ
C      IN THE BUFFERS IN ONE CALL AND THEN RETURN.
C
C      FLOW:
C
C      1. SET NUMBER OF BLOCKS TO SEARCH FOR.  SET THE RECALL FLAG
C         = .FALSE.
C
C      2. DO A GLOBAL SWEEP OF ALL BUFFERS LOOKING FOR THE NECESSARY
C         WORD ADDRESSES IN ANY ORDER.
C
C      3. IF WE FOUND THEM ALL EXIT.
C
C      4. SEARCH *WABP* ARRAY FOR BLOCK NUMBER TO MATCH FIRST.
C
C      5. IF NO MATCH CALL *WGETB* AND THEN RETURN.
C
C      6. WE HAVE FOUND THE FIRST BLOCK. IF IT IS THE LAST PHYSICAL
C         ONE IN THE BUFFER POOL THAN SKIP TO 10.
C
C      7. LOOK AT NEXT BLOCK.  IF IT DOESNT MATCH THE NEXT DESIRED
C         WORD ADDRESS SKIP TO 10.
C
C      8. DECREMENT BLOCK COUNTER.
C
C      9. IF ALL BLOCKS NOT SEARCHED LOOP TO 7.
C
C     10. WE HAVE FOUND AT LEAST ONE BLOCK. IF THE BLOCK OR BLOCKS
C         WERE THE LAST PHYSICAL ONES IN THE BUFFER POOL SET THE
C         RECALL FLAG = .TRUE. AND RETURN.
C
C     11. WE ARE SHORT BLOCKS AND THERE IS ROOM AT THE END OF THE
C         BUFFER POOL.  IF WE ARE READING RECOMPUTE THE ENDING
C         BLOCK NUMBER BASED ON THE HISTORY OF SEQUENTIAL READING
C         DONE UP TO THIS POINT.  IF THERE IS MORE THAN ONE BLOCK
C         INVOLVED OR IF IT IS ONE BLOCK AND IT IS IN WRITE MODE
C         *WFLUSH* IS CALLED TO FLUSH TO DISK.
C
C     12. PUT THE WORD ADDRESSES OF THE BUFFERS WE NEED IN THE
C         *WABP* ARRAY, MAKING SURE THERE ARE NO DUPLICATE ENTRIES.
C
C     13. IF READING FROM DISK SKIP TO 16.
C
C     14. IF WRITING AND THE LAST BLOCK ADDRESS PUT INTO THE *WABP*
C         ARRAY IS THE LAST BLOCK TO TRANSFER AND IF THIS LAST
C         BLOCK IS NOT BEYOND EOI AND IF THIS LAST BLOCK IS A
C         PARTIAL ONE THEN WE MUST READ IT IN FROM DISK TO MERGE
C         THE OLD AND THE NEW DATA TOGETHER.
C
C     15. SKIP TO 17.
C
C     16. READ IN THE BLOCKS TO THE HIGH BUFFERS.
C
C     17. RETURN.

CDIR$  EJECT
       IMPLICIT INTEGER (A-Z)

       INCLUDE "wavars.fh"

       LOGICAL RECALL

       WABP(BUFFERS+1) = 0
       RECALL = .FALSE.
       BLOCKO = MIN0(BUFFERS,BLKCNT)
       K      = BLOCKO

C      Search for the first BLOCKO file pages in the buffer cache.
C      If all are present, then we return.   An inital segment of the
C      current request is completely buffer resident in that case.

       WADD1   = BSIZE * (FIRST - 1) + 1
       BK      = K

       DO  5 I = 1,K
           IF (IXMM@(WABP,BUFFERS,MSK,WADD1).NE.0) BK = BK - 1
           WADD1   = WADD1 + BSIZE
5      CONTINUE
       IF(BK.EQ.0) THEN
           FET(4,INDEX) = FET(4,INDEX) + SHIFTL(1,32)
           RETURN
       ENDIF

C      If the large initial segment was not completely buffer resident, then
C      rescan to see if at least the first page is buffer resident.

       WADD1  = BSIZE * (FIRST - 1) + 1
       IBLK   = FIRST
       HIGH   = 0
       LOW    = BUFFERS + 1
       BEGIN  = .FALSE.
       I      = IXMM@(WABP,BUFFERS,MSK,WADD1)

       IF(I.EQ.0) THEN

C          The first word address is not buffer resident.   Call WGETB to
C          assign a buffer to this page (and those pages which follow?).

           CALL WGETB(0)
           FET(3,INDEX) = FET(3,INDEX) + SHIFTL(1,32)
           RETURN
       ENDIF

       K    = K - 1
       ISA  = I + 1
       IBLK = FIRST+1

C      Find the largest initial contiguous segment which is buffer resident.

       IF (I.NE.BUFFERS) THEN
           DO 10 J = I+1,BUFFERS+1
               WADD1   = WADD1 + BSIZE
               IF ((WABP(J).AND.MSK).EQ.WADD1) THEN
                   K    = K - 1
                   ISA  = J + 1
                   IBLK = IBLK+1
               ELSE
                   GO TO 20
               ENDIF
10         CONTINUE
20         CONTINUE
       ENDIF

C      Here we have found some buffer matches but not all.

       FET(5,INDEX) = FET(5,INDEX) + SHIFTL(1,32)       ! partial hit count
       IF (I.EQ.BUFFERS .OR. J.EQ.(BUFFERS+1)) THEN
           RECALL = .TRUE.                              ! at end of buffer
           RETURN
       ENDIF

C      Check if there's room at the end for more file pages for this request.

       SHORT = K
       END   = ISA + SHORT - 1
       LIM   = ISA + (EOIB-IBLK)
       IF (READ) THEN
           LSTSEQ        = FET(12,INDEX).AND.MASK(128-32)
           LSTSEQ        = (LSTSEQ + SHIFTR(FET(11,INDEX),32)) / 2
           FET(12,INDEX) = FET(12,INDEX).AND.MASK(32) + LSTSEQ
           END           = MAX0(END,ISA + LSTSEQ - 1)
           END           = MIN0(END,LIM)
       ENDIF
       IF (END.GT.BUFFERS) END = BUFFERS

C      Check for dirty blocks between ISA and END inclusive.   If there are
C      any, then flush them to disk.  The signbit indicates dirty buffers.

       IF (END.GT.(ISA+1)) THEN
           IX = IXMM@(WABP(ISA),END-ISA+1,SIGNBIT,SIGNBIT)
           IF(IX.NE.0) CALL WFLUSH(ISA,END)
       ELSE
           IF((WABP(ISA).LT.0).OR.(WABP(END).LT.0))CALL WFLUSH(ISA,END)
       ENDIF

C      Check for trailing pages which are buffer resident in another part of
C      the buffer.   Invalidate these pages.  We will assign new blocks
C      to these file pages so that the initial segment is contiguous in
C      the buffer.

       WADD1   = (WABP(ISA-1).AND.MSK)
       IF ((EOIB.LT.IBLK).AND.WRITE) THEN
           DO 30 I = ISA,END
               WADD1   = WADD1 + BSIZE
               WABP(I) = WADD1
30         CONTINUE
       ELSE
           DO 40 I = ISA,END
               WADD1   = WADD1 + BSIZE
               J       = IXMM@(WABP,BUFFERS,MSK,WADD1)
               IF (J.NE.0) THEN
                   IF (WABP(J).LT.0) THEN
                       CALL WFLUSH(1,BUFFERS)
                   ENDIF
                   WABP(J) = 2                          ! mark page as unused
               ENDIF
               WABP(I) = WADD1
40         CONTINUE
       ENDIF

       IF (WRITE) THEN
           IF (BLKCNT.LE.BUFFERS .AND. PARTEND .AND. WADD1.EQ.LBWA .AND.
     +         LAST.LE.EOIB) THEN

C              We are processing the last block containing the tail of the
C              requested range of word addresses.  If the request ends between
C              block boundaries, the block must be pre-read from disk.
C
C              But first flush the rest of the buffer to disk if the block we
C              are about to read would otherwise be past the end of file.

               IF (LAST.GT.EOIBW) CALL WFLUSH(1,BUFFERS)

               II    = BSIZE * (END - 1) + 1
               FET(6,INDEX) = FET(6,INDEX) + 1      ! cnt of READWA calls
               FET(13,INDEX)= FET(13,INDEX) + BSIZE ! wrd cnt of phys xfers
               IT1          = IRTC()
               CALL READWA(INDEX,WAA(II),WADD1,BSIZE,0)
               FET(10,INDEX)= FET(10,INDEX)+(IRTC()-IT1) ! I/O time
           ENDIF
       ELSE

C          Read in the following sequence of file blocks.   But first flush
C          dirty buffers if necessary to extend the disk image to the full
C          file size.

           WADDR = WABP(ISA)
           II    = BSIZE * (ISA - 1) + 1
           ICDK  = BSIZE * (END - ISA + 1)
           IF (((WADDR+ICDK-1)/BSIZE).GT.EOIBW) CALL WFLUSH(1,BUFFERS)
           FET(6,INDEX) = FET(6,INDEX) + 1          ! cnt of READWA calls
           IT1          = IRTC()
           FET(13,INDEX)= FET(13,INDEX) + ICDK      ! wrd cnt of phys xfers

           CALL READWA(INDEX,WAA(II),WADDR,ICDK,0)
           FET(10,INDEX) = FET(10,INDEX) + (IRTC() - IT1)

       ENDIF

       RETURN
       END

ext: f

max_stars_repo_name: sharugupta/OpenUH

max_issues_repo_name: sharugupta/OpenUH

max_forks_repo_name: sharugupta/OpenUH


text:
IT Express is the computing help desk for students, faculty and staff at UC Davis. They are available via phone (530) 7544357 and via email @ ithelp@ucdavis.edu.

A list of supported services can be found http://itexpress.ucdavis.edu/support/ here.

The IT Express Knowledge base or xbase is located at http://xbase.ucdavis.edu http://xbase.ucdavis.edu.

If you are looking for computer repair on campus, seek out IT Professional Services.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
** Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
** See https://llvm.org/LICENSE.txt for license information.
** SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

*   ABS,MIN,MAX constant folding(mostly) in the expander.

	parameter(n = 83)
	integer rslts(n), expect(n)
	real  rrslts(n), rexpect(n)
	equivalence (rslts, rrslts), (expect, rexpect)

        integer  i, j, k, l
        real     r, s, t, u
        double precision  d, e, f, g

	data (expect(i), i = 1, 5) /2, 34, 0, 34, 2 /
        i = ABS( 2)
        rslts(1) = i
        j = ABS( -34)
        rslts(2) = j
        k = ABS(0)
        rslts(3) = k
        l = ABS(j)
        rslts(4) = l
        l = ABS(i)
        rslts(5) = l

	data (rexpect(i), i = 6, 9) /2.45, 34.6, 0.0, 34.6/
        r = ABS(2.45)
        rrslts(6) = r
        s = ABS(-34.6)
        rrslts(7) = s
        t = ABS(0.0)
        rrslts(8) = t
        u = ABS(s)
        rrslts(9) = u

	data (rexpect(i), i = 10, 13) /2.45, 34.6, 0.0, 34.6/
        d = ABS(2.45D0)
        rrslts(10) = d
        e = ABS(-34.6D0)
        rrslts(11) = e
        f = ABS(0.0D0)
        rrslts(12) = f
        g = ABS(e)
        rrslts(13) = g

        data (expect(i), i = 14, 24) /23,23,-23,-23,23,34, 23,
     - 23, 0 ,34, 34/
        i = MAX( 2, 23)
        rslts(14) = i
        i = MAX( 23, 2)
        rslts(15) = i
        j = MAX( -34, -23)
        rslts(16) = j
        j = MAX( -23, -34)
        rslts(17) = j
        j = MAX( 23, -34)
        rslts(18) = j
        j = MAX( -23, 34)
        rslts(19) = j
        k = MAX(0, 23)
        rslts(20) = k
        k = MAX(23, 0)
        rslts(21) = k
        k = MAX(0,0)
        rslts(22) = k
        l = MAX(j,k)
        rslts(23) = l
        l = MAX(k,l)
        rslts(24) = l

	data (rexpect(i), i = 25, 36) /23.2,23.2,-23.4,-23.4, 34.6,
     - 23.4, 0.0, 34.5, 0.0, 23.4, 23.4, 23.4 /
        r = MAX(2.45,23.2)
        rrslts(25) = r
        r = MAX(23.2,2.45)
        rrslts(26) = r
        s = MAX(-34.6,-23.4)
        rrslts(27) = s
        s = MAX(-23.4,-34.6)
        rrslts(28) = s
        s = MAX(-23.4,34.6)
        rrslts(29) = s
        s = MAX(23.4,-34.6)
        rrslts(30) = s
        t = MAX(0.0, -34.5)
        rrslts(31) = t
        t = MAX(0.0, 34.5)
        rrslts(32) = t
        t = MAX(0.0, 0.0)
        rrslts(33) = t
        u = MAX(s,t)
        rrslts(34) = u
        u = MAX(t,s)
        rrslts(35) = u
        u = MAX(s,s)
        rrslts(36) = u

        data (rexpect(i),i=37,48)/34.3D0,34.3D0,-23.6D0,-34.6D0,
     - 34.6D0,57.9D0,34.8D0,34.8D0,0.0D0,57.9D0, 57.9D0, 57.9D0/
        d = MAX(2.45D0,34.3D0)
        rrslts(37) = d
        d = MAX(34.3D0, 2.45D0)
        rrslts(38) = d
        e = MAX(-34.6D0,-23.6D0)
        rrslts(39) = e
        e = MAX(-34.6D0, -57.9D0)
        rrslts(40) = e
        e = MAX(34.6D0, -57.9D0)
        rrslts(41) = e
        e = MAX(-34.6D0, 57.9D0)
        rrslts(42) = e
        f = MAX(0.0D0,34.8D0)
        rrslts(43) = f
        f = MAX(34.8D0,0.0D0)
        rrslts(44) = f
        f = MAX(0.0D0,0.0D0)
        rrslts(45) = f
        g = MAX(e,f)
        rrslts(46) = g
        g = MAX(g,e)
        rrslts(47) = g
        g = MAX(g,g)
        rrslts(48) = g

	data (expect(i),i= 49,59)/2,2,-34,-34,-34,-23,0,0,0,-23,-23/
        i = MIN( 2, 23)
        rslts(49) = i
        i = MIN( 23, 2)
        rslts(50) = i
        j = MIN( -34, -23)
        rslts(51) = j
        j = MIN( -23, -34)
        rslts(52) = j
        j = MIN( 23, -34)
        rslts(53) = j
        j = MIN( -23, 34)
        rslts(54) = j
        k = MIN(0, 23)
        rslts(55) = k
        k = MIN(23, 0)
        rslts(56) = k
        k = MIN(0,0)
        rslts(57) = k
        l = MIN(j,k)
        rslts(58) = l
        l = MIN(k,l)
        rslts(59) = l

	data (rexpect(i), i = 60, 71) /2.45,2.45,-34.6,-34.6,-23.4,-34.6,
     - -34.5, 0.0,0.0,-34.6,-34.6,-34.6/
        r = MIN(2.45,23.2)
        rrslts(60) = r
        r = MIN(23.2,2.45)
        rrslts(61) = r
        s = MIN(-34.6,-23.4)
        rrslts(62) = s
        s = MIN(-23.4,-34.6)
        rrslts(63) = s
        s = MIN(-23.4,34.6)
        rrslts(64) = s
        s = MIN(23.4,-34.6)
        rrslts(65) = s
        t = MIN(0.0, -34.5)
        rrslts(66) = t
        t = MIN(0.0, 34.5)
        rrslts(67) = t
        t = MIN(0.0, 0.0)
        rrslts(68) = t
        u = MIN(s,t)
        rrslts(69) = u
        u = MIN(t,s)
        rrslts(70) = u
        u = MIN(s,s)
        rrslts(71) = u

	data(rexpect(i),i = 72, 83)/2.45D0,2.45D0,-34.6D0,-57.9D0,
     - -57.9D0,-34.6D0,0.0D0,0.0D0,0.0D0,-34.6D0,-34.6D0,-34.6D0/
        d = MIN(2.45D0,34.3D0)
        rrslts(72) = d
        d = MIN(34.3D0, 2.45D0)
        rrslts(73) = d
        e = MIN(-34.6D0,-23.6D0)
        rrslts(74) = e
        e = MIN(-34.6D0, -57.9D0)
        rrslts(75) = e
        e = MIN(34.6D0, -57.9D0)
        rrslts(76) = e
        e = MIN(-34.6D0, 57.9D0)
        rrslts(77) = e
        f = MIN(0.0D0,34.8D0)
        rrslts(78) = f
        f = MIN(34.8D0,0.0D0)
        rrslts(79) = f
        f = MIN(0.0D0,0.0D0)
        rrslts(80) = f
        g = MIN(e,f)
        rrslts(81) = g
        g = MIN(g,e)
        rrslts(82) = g
        g = MIN(g,g)
        rrslts(83) = g

	call check(rslts, expect, n)

        end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
! fs5_def.for

	parameter	w_fs5_mode_year  =31
	parameter	w_fs5_mode_month =32
	parameter	w_fs5_mode_dom   =33


text:
subroutine	mcsigma(chiscale,w,pf,pf0,obs,lambda_obs,mobs,lsfarr,e,sp,lchi,cov)

!
!	Monte-Carlo error estimates of the error bars.
!
!	We search for solutions after introducing Gaussian errors
!	in the observations and derive the standard deviation for
!	each parameter.
!

use random
use share, only: dp,ndim,nov,indv,nlambda1,inter,algor,mcruns,mlsf,nlsf

implicit none

!in/out
real(dp), intent(in)	   :: chiscale		!chi**2/sum(w_i(f_i-obs_i))^2)
real(dp), intent(in)       :: w(nlambda1)       ! weights
real(dp), intent(in)       :: pf(ndim)          ! vector of fixed parameters
real(dp), intent(in)       :: pf0(ndim)         ! pars read from pfile (in physical units)
real(dp), intent(inout)    :: obs(nlambda1)     ! vector of observations
real(dp), intent(in)       :: lambda_obs(nlambda1)  ! vector of wavelength for observations
real(dp), intent(in)       :: mobs                ! mean or median of obs array
real(dp), intent(in)       :: lsfarr(mlsf,nlsf)    ! lsfarray
real(dp), intent(in)  	   :: e(nlambda1)		!errors in the observations
real(dp), intent(out)	   :: sp(ndim)		!vector of uncertainties
real(dp), intent(out)  	   :: lchi		!log10(chi**2/(npix-nov+1))
real(dp), intent(out)      :: cov(nov,nov)	!covariance matrix

!locals
real(dp)                   :: pf2(ndim)         !vector of fixed parameters that will get modified
						!in individual calls to getmin
integer		  	   :: i,j,k		!counters
real(dp)		   :: p(nov)		!vector of var. pars. (1 run)
real(dp)		   :: pm(nov)		!average values
real(dp)		   :: pt(mcruns,nov) !tmp storage for the results
real(dp) 	  	   :: hobs(nlambda1) !holder for the original obs
real(dp)		   :: flux(nlambda1) !interp. flux vector

!initialize
pm(:)=0.0_dp
cov(:,:)=0.0_dp
sp(:)=0.0_dp

!keep the original data; the obs array will be modified here
hobs=obs

!evaluate lchi

call flx(pf,lambda_obs,e,mobs,lsfarr,flux)
lchi=sum(w*(obs(:)-flux(:))**2)
lchi=log10(lchi*chiscale/(nlambda1-nov+1))

if (lchi >= 4.0d0) then
!don't even try with lost causes. Very poorly matched spectra
!are simply assigned sp(*)=-1.00

  	do j=1,nov
		sp(indv(j))=-1.00_dp
  	enddo

else

	!initialize the random number generator
	call random_seed()

	!repeat mcruns times
	do j=1,mcruns

		!introduce Gaussian noise
		do i=1,nlambda1
			obs(i)=random_normal()*e(i)+hobs(i)
		enddo

		!initialize
		p(:)=0.5_dp
		!call random_number(p)

		!solve
		pf2=pf
		call getmin(algor,w,pf2,pf0,obs,lambda_obs,e,mobs,lsfarr,p)

		!store info
		pt(j,:)=p(:)  	!store results to calculate covariance later

		!write(*,*)'mcrun number/sol =',j,p

	enddo

	!compute mean
	do j=1,nov
		pm(j)=sum(pt(:,j))/mcruns
	enddo
	!and covariance
	do j=1,nov
	  do k=1,nov
 	   cov(k,j)=sum( (pt(:,k)-pm(k))*(pt(:,j)-pm(j)) )/(mcruns-1._dp)
	  enddo
	enddo

	!use diagonal elements to get std. deviation
  	do j=1,nov
		if (cov(j,j) >=0.0) then
			sp(indv(j))=sqrt(cov(j,j))
		else
			sp(indv(j))=-1._dp
		endif
  	enddo

endif

!return the original data to obs
obs=hobs

end subroutine	mcsigma

max_stars_repo_name: dnidever/apogee

max_issues_repo_name: dnidever/apogee

max_forks_repo_name: dnidever/apogee


text:
Now home to Barnyard Theatre, this historic structure was built by German immigrant Gottfried Schmeiser in the 1860s. The story of the barn, as written by Carla Huntone of the current Schmeiser Ranch residents and property ownersis as follows:

In the mid1800s, Gottfried Schmeiser, a young carpenter, left Germany for the New World. After obtaining a land grant signed by Ulysses S. Grant, he headed west to settle his land. By 1868, he had built a lovely home, planted fields of wheat and an almond orchard, and married his neighbor, Sophia Oeste, for which Oeste Drive is named. (See also Theodore Schmeiser.)

Gottfried and Sophia Schmeiser raised six children in the house that still stands on the property. Only one of the six children wanted to continue the family farm, and that child was a daughter! As unusual as it was at that time for a single woman to run a farm, Elmine Schmeiser did a bangup job. Some Yolo County oldtimers still remember Elmine, who lived well into her eighties and ran the 160acre Schmeiser Ranch.

The ranch was next taken care of by neighbor Doc Russell, whose ranch house is still visible just south of Schmeiser Ranch. When Doc Russell passed on, he willed the land to UCD for agricultural research. The fields surrounding the homestead are owned by UCD, and leased to wellknown Yolo County farmer Jack Meek.

For some time, the beautiful house and barns were rented out to students, and the homestead steadily declined. The current barn is the new barn, built in 1893. The old barn with its attached water cistern burned to the ground while students still rented the land. Fortunately, the home was saved. Soon after, Tosi Marceline and Phil Kitchen rented the homestead and started a family and a business there, the Three Palms Nursery, which is still in operation at Russell Boulevard and County Road 95A.

For 20 years, Phil and Tosi loved the homestead and cared for it, making historically accurate renovations and maintaining the lovely home and barn. In 1996, the university put the house and nine surrounding acres up for sale. Charles and Carla Hunt (who, like the Schmeisers before them, have six children), bought the Schmeiser Ranch and continue to treasure it. The current barn is operational, filled with horses, cows, pigs, goats, sheep, barn cats, chickens and 4H kids.