osse/EnKF/spofa.F

  
  
       subroutine spofa(a,lda,n,info)
      integer lda,n,info
      real a(lda,1)
c
c     spofa factors a real symmetric positive definite matrix.
c
c     spofa is usually called by spoco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for spoco) = (1 + 18/n)*(time for spofa) .
c
c     on entry
c
c        a       real(lda, n)
c                the symmetric matrix to be factored.  only the
c                diagonal and upper triangle are used.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix  r  so that  a = trans(r)*r
c                where  trans(r)  is the transpose.
c                the strict lower triangle is unaltered.
c                if  info .ne. 0 , the factorization is not complete.
c
c        info    integer
c                = 0  for normal return.
c                = k  signals an error condition.  the leading minor
c                     of order  k  is not positive definite.
c
c     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas sdot
c     fortran sqrt
c
c     internal variables
c
      real sdot,t
      real s
      integer j,jm1,k
c     begin block with ...exits to 40
c
c
         do 30 j = 1, n
            info = j
            s = 0.00
            jm1 = j - 1
            if (jm1 .lt. 1) go to 20
            do 10 k = 1, jm1
               t = a(k,j) - sdot(k-1,a(1,k),1,a(1,j),1)
               t = t/a(k,k)
               a(k,j) = t
               s = s + t*t
   10       continue
   20       continue
            s = a(j,j) - s
c     ......exit
            if (s .le. 0.00) go to 40
            a(j,j) = sqrt(s)
   30    continue
         info = 0
   40 continue
      return
      end
1	afe	1.1
2
3			subroutine spofa(a,lda,n,info)
4			integer lda,n,info
5			real a(lda,1)
6			c
7			c spofa factors a real symmetric positive definite matrix.
8			c
9			c spofa is usually called by spoco, but it can be called
10			c directly with a saving in time if rcond is not needed.
11			c (time for spoco) = (1 + 18/n)*(time for spofa) .
12			c
13			c on entry
14			c
15			c a real(lda, n)
16			c the symmetric matrix to be factored. only the
17			c diagonal and upper triangle are used.
18			c
19			c lda integer
20			c the leading dimension of the array a .
21			c
22			c n integer
23			c the order of the matrix a .
24			c
25			c on return
26			c
27			c a an upper triangular matrix r so that a = trans(r)*r
28			c where trans(r) is the transpose.
29			c the strict lower triangle is unaltered.
30			c if info .ne. 0 , the factorization is not complete.
31			c
32			c info integer
33			c = 0 for normal return.
34			c = k signals an error condition. the leading minor
35			c of order k is not positive definite.
36			c
37			c linpack. this version dated 08/14/78 .
38			c cleve moler, university of new mexico, argonne national lab.
39			c
40			c subroutines and functions
41			c
42			c blas sdot
43			c fortran sqrt
44			c
45			c internal variables
46			c
47			real sdot,t
48			real s
49			integer j,jm1,k
50			c begin block with ...exits to 40
51			c
52			c
53			do 30 j = 1, n
54			info = j
55			s = 0.00
56			jm1 = j - 1
57			if (jm1 .lt. 1) go to 20
58			do 10 k = 1, jm1
59			t = a(k,j) - sdot(k-1,a(1,k),1,a(1,j),1)
60			t = t/a(k,k)
61			a(k,j) = t
62			s = s + t*t
63			10 continue
64			20 continue
65			s = a(j,j) - s
66			c ......exit
67			if (s .le. 0.00) go to 40
68			a(j,j) = sqrt(s)
69			30 continue
70			info = 0
71			40 continue
72			return
73			end