CMAX: [Version 1.0 (6/22/93)]
CMAX: ----------------------------------------------------
CMAX: Transformation of VV2 from vv.f
CMAX: ----------------------------------------------------
C----------------------------------------------------------------------
C
C  CHOLESKY DECOMPOSITION OF POSITIVE DEFINITE MATRIX OF ORDER NR
C  WITH SOLUTION OF NQ COLUMN VECTORS
C
C  A(NS,NR*(NR+1)/2) = POSITIVE DEFINITE MATRIX STACK (TRIANGULAR FORMS)
C  B(NS,NR,NQ)       = RHS MATRIX STACK (INPUT)
C  X(NS,NR,NQ)       = SOLUTION MATRIX STACK (OUTPUT)
C  DINV,T            = WORKING ARRAYS
C  MARK              = ON OUTPUT CONTAINS THE NUMBER OF DIAGONAL
C                      ELEMENTS SUCCESSFULY DECOMPOSED IN EACH MATRIX
C  NIN               = NUMBER OF MATRICES IN INPUT STACK
C  NMAX              = LARGEST MATRIX IN STACK
C  NS,NR             = DIMENSIONS OF MATRIX STACKS
C  NQ                = WIDTH OF FORCING MATRIXES
C
C  TRIANGULAR REPRESENTATION LISTS ELEMENTS TOWARDS THE DIAGONAL
C
C  Program history:
C     ?/??/??   Algorithm written            H. Carus
C     3/31/84   Original code                J. Woollen
C    11/14/90   Recoded for Cray YMP         M. Seablom
C     2/13/91   Tolerance changed to 1.e-3   M. Seablom
C     3/05/93   Removal of matrix stacks     M. Seablom
C
C----------------------------------------------------------------------
C                                                                     !
C NASA Goddard Space Flight Center       Data Assimilation Office     !
C                                                                     !
C Goddard Earth Observing System                                      !
C Data Analysis System - Version 1.5                                  !
C October, 1993                                                       !
C                                                                     !
c---------------------------------------------------------------------c
      SUBROUTINE vv2(a,b,dmark,nh,nin,nq,nr,ml,ns,nv,x,rdummy4)
      IMPLICIT REAL (a - h, o - z)
      REAL a (ns, ml), b (ns, nr, nv, nh, nq), x (ns, nr, nv, nh, nq)
      REAL t, vars, dmark
!hpf$ distribute a(block,block)
!hpf$ distribute b(block,block,block,block,block)
!hpf$ distribute x(block,block,block,block,block)
      DATA icounter / 0 /
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C  Initialize arrays.
      ix(i,j) = i * (i - 1) / 2 + j
      IF (nin .NE. 1) THEN
        PRINT *, 'nin is not equal to 1 ',nin
        STOP
      ENDIF
      icounter = icounter + 1
      t = 0.
      vars = 0.
      n = nr
C  Decompose the matrices.
      DO 7004 i = 1 , n
        DO 7002 j = 1 , i
          t = a(1,ix(i,j))
          vars = a(1,ix(i,j))
          jm1 = j - 1
          IF (j .NE. 1) THEN
            DO 7000 k = 1 , jm1
              t = t - a(1,ix(i,k)) * a(1,ix(j,k))
 7000       CONTINUE
          ENDIF
          IF (i .GT. j) THEN
            a(1,ix(i,j)) = t * a(1,ix(j,j))
          ENDIF
          IF (i .EQ. j) THEN
            IF ((t.GT.(1.E-3 * vars)).AND.(dmark .EQ. float(i - 1))
     &            ) THEN
              a(1,ix(i,i)) = 1. / sqrt(t)
              dmark = float(i)
            ELSE
              a(1,ix(i,i)) = 0.0
            ENDIF
          ENDIF
 7002   CONTINUE
 7004 CONTINUE
      n = j - 1
C  Solve upper and lower triangular systems for all right hand sides.
      DO 7014 nn = 1 , nv * nh * nq
C  Forward substitution.
        DO 7008 i = 1 , n
          t = b(1,i,nn,1,1)
          IF (i .NE. 1) THEN
            DO 7006 j = 1 , (i - 1)
              t = t - a(1,ix(i,j)) * x(1,j,nn,1,1)
 7006       CONTINUE
          ENDIF
          x(1,i,nn,1,1) = t * a(1,ix(i,i))
 7008   CONTINUE
C  Backward substitution.
        DO 7012 ii = 1 , n
          i = n + 1 - ii
          t = x(1,i,nn,1,1)
          IF (i .NE. n) THEN
            DO 7010 j = i + 1 , n
              t = t - a(1,ix(j,i)) * x(1,j,nn,1,1)
 7010       CONTINUE
          ENDIF
          x(1,i,nn,1,1) = t * a(1,ix(i,i))
 7012   CONTINUE
C  End of solution for each forcing.
 7014 CONTINUE
      RETURN
      END
