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