osse/EnKF/sposl.F

   
  
  
  
        subroutine sposl(a,lda,n,b)
      integer lda,n
      real a(lda,1),b(1)
c
c     sposl solves the real symmetric positive definite system
c     a * x = b
c     using the factors computed by spoco or spofa.
c
c     on entry
c
c        a       real(lda, n)
c                the output from spoco or spofa.
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        b       real(n)
c                the right hand side vector.
c
c     on return
c
c        b       the solution vector  x .
c
c     error condition
c
c        a division by zero will occur if the input factor contains
c        a zero on the diagonal.  technically this indicates
c        singularity but it is usually caused by improper subroutine
c        arguments.  it will not occur if the subroutines are called
c        correctly and  info .eq. 0 .
c
c     to compute  inverse(a) * c  where  c  is a matrix
c     with  p  columns
c           call spoco(a,lda,n,rcond,z,info)
c           if (rcond is too small .or. info .ne. 0) go to ...
c           do 10 j = 1, p
c              call sposl(a,lda,n,c(1,j))
c        10 continue
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 saxpy,sdot
c
c     internal variables
c
      real sdot,t
      integer k,kb
c
c     solve trans(r)*y = b
c
      do 10 k = 1, n
         t = sdot(k-1,a(1,k),1,b(1),1)
         b(k) = (b(k) - t)/a(k,k)
   10 continue
c
c     solve r*x = y
c
      do 20 kb = 1, n
         k = n + 1 - kb
         b(k) = b(k)/a(k,k)
         t = -b(k)
         call saxpy(k-1,t,a(1,k),1,b(1),1)
   20 continue
      return
      end
1	afe	1.1
2
3
4
5			subroutine sposl(a,lda,n,b)
6			integer lda,n
7			real a(lda,1),b(1)
8			c
9			c sposl solves the real symmetric positive definite system
10			c a * x = b
11			c using the factors computed by spoco or spofa.
12			c
13			c on entry
14			c
15			c a real(lda, n)
16			c the output from spoco or spofa.
17			c
18			c lda integer
19			c the leading dimension of the array a .
20			c
21			c n integer
22			c the order of the matrix a .
23			c
24			c b real(n)
25			c the right hand side vector.
26			c
27			c on return
28			c
29			c b the solution vector x .
30			c
31			c error condition
32			c
33			c a division by zero will occur if the input factor contains
34			c a zero on the diagonal. technically this indicates
35			c singularity but it is usually caused by improper subroutine
36			c arguments. it will not occur if the subroutines are called
37			c correctly and info .eq. 0 .
38			c
39			c to compute inverse(a) * c where c is a matrix
40			c with p columns
41			c call spoco(a,lda,n,rcond,z,info)
42			c if (rcond is too small .or. info .ne. 0) go to ...
43			c do 10 j = 1, p
44			c call sposl(a,lda,n,c(1,j))
45			c 10 continue
46			c
47			c linpack. this version dated 08/14/78 .
48			c cleve moler, university of new mexico, argonne national lab.
49			c
50			c subroutines and functions
51			c
52			c blas saxpy,sdot
53			c
54			c internal variables
55			c
56			real sdot,t
57			integer k,kb
58			c
59			c solve trans(r)*y = b
60			c
61			do 10 k = 1, n
62			t = sdot(k-1,a(1,k),1,b(1),1)
63			b(k) = (b(k) - t)/a(k,k)
64			10 continue
65			c
66			c solve r*x = y
67			c
68			do 20 kb = 1, n
69			k = n + 1 - kb
70			b(k) = b(k)/a(k,k)
71			t = -b(k)
72			call saxpy(k-1,t,a(1,k),1,b(1),1)
73			20 continue
74			return
75			end