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
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