The R Project SVN R

Rev

Rev 3388 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 3388 Rev 12778
Line 1... Line 1...
1 
      subroutine dpbfa(abd,lda,n,m,info)
 1 
      subroutine dpbfa(abd,lda,n,m,info)
2 
 
 2 
 
3 
      implicit none
 - 
 
4 
      integer lda,n,m,info
 3 
      integer lda,n,m,info
5 
      double precision abd(lda,n)
 4 
      double precision abd(lda,n)
6 
c
 5 
c
7 
c     dpbfa factors a double precision symmetric positive definite
 6 
c     dpbfa factors a double precision symmetric positive definite
8 
c     matrix stored in band form.
 7 
c     matrix stored in band form.
9 
c
 8 
c
10 
c     dpbfa is usually called by dpbco, but it can be called
 9 
c     dpbfa is usually called by dpbco, but it can be called
11 
c     directly with a saving in time if	 rcond	is not needed.
 10 
c     directly with a saving in time if  rcond  is not needed.
12 
c
 11 
c
13 
c     on entry
 12 
c     on entry
14 
c
 13 
c
15 
c	 abd	 double precision(lda, n)
 14 
c        abd     double precision(lda, n)
16 
c		 the matrix to be factored.  the columns of the upper
 15 
c                the matrix to be factored.  the columns of the upper
17 
c		 triangle are stored in the columns of abd and the
 16 
c                triangle are stored in the columns of abd and the
18 
c		 diagonals of the upper triangle are stored in the
 17 
c                diagonals of the upper triangle are stored in the
19 
c		 rows of abd .	see the comments below for details.
 18 
c                rows of abd .  see the comments below for details.
20 
c
 19 
c
21 
c	 lda	 integer
 20 
c        lda     integer
22 
c		 the leading dimension of the array  abd .
 21 
c                the leading dimension of the array  abd .
23 
c		 lda must be .ge. m + 1 .
 22 
c                lda must be .ge. m + 1 .
24 
c
 23 
c
25 
c	 n	 integer
 24 
c        n       integer
26 
c		 the order of the matrix  a .
 25 
c                the order of the matrix  a .
27 
c
 26 
c
28 
c	 m	 integer
 27 
c        m       integer
29 
c		 the number of diagonals above the main diagonal.
 28 
c                the number of diagonals above the main diagonal.
30 
c		 0 .le. m .lt. n .
 29 
c                0 .le. m .lt. n .
31 
c
 30 
c
32 
c     on return
 31 
c     on return
33 
c
 32 
c
34 
c	 abd	 an upper triangular matrix  r , stored in band
 33 
c        abd     an upper triangular matrix  r , stored in band
35 
c		 form, so that	a = trans(r)*r .
 34 
c                form, so that  a = trans(r)*r .
36 
c
 35 
c
37 
c	 info	 integer
 36 
c        info    integer
38 
c		 = 0  for normal return.
 37 
c                = 0  for normal return.
39 
c		 = k  if the leading minor of order  k	is not
 38 
c                = k  if the leading minor of order  k  is not
40 
c		      positive definite.
 39 
c                     positive definite.
41 
c
 40 
c
42 
c     band storage
 41 
c     band storage
43 
c
 42 
c
44 
c	    if	a  is a symmetric positive definite band matrix,
 43 
c           if  a  is a symmetric positive definite band matrix,
45 
c	    the following program segment will set up the input.
 44 
c           the following program segment will set up the input.
46 
c
 45 
c
47 
c		    m = (band width above diagonal)
 46 
c                   m = (band width above diagonal)
48 
c		    do 20 j = 1, n
 47 
c                   do 20 j = 1, n
49 
c		       i1 = max0(1, j-m)
 48 
c                      i1 = max0(1, j-m)
50 
c		       do 10 i = i1, j
 49 
c                      do 10 i = i1, j
51 
c			  k = i-j+m+1
 50 
c                         k = i-j+m+1
52 
c			  abd(k,j) = a(i,j)
 51 
c                         abd(k,j) = a(i,j)
53 
c		 10    continue
 52 
c                10    continue
54 
c		 20 continue
 53 
c                20 continue
55 
c
 54 
c
56 
c     linpack.	this version dated 08/14/78 .
 55 
c     linpack.  this version dated 08/14/78 .
57 
c     cleve moler, university of new mexico, argonne national lab.
 56 
c     cleve moler, university of new mexico, argonne national lab.
58 
c
 57 
c
59 
c     subroutines and functions
 58 
c     subroutines and functions
60 
c
 59 
c
61 
c     blas ddot
 60 
c     blas ddot
Line 67... Line 66...
67 
      double precision s
 66 
      double precision s
68 
      integer ik,j,jk,k,mu
 67 
      integer ik,j,jk,k,mu
69 
c     begin block with ...exits to 40
 68 
c     begin block with ...exits to 40
70 
c
 69 
c
71 
c
 70 
c
72 
	 do 30 j = 1, n
 71 
         do 30 j = 1, n
73 
	    info = j
 72 
            info = j
74 
	    s = 0.0e0
 73 
            s = 0.0e0
75 
	    ik = m + 1
 74 
            ik = m + 1
76 
	    jk = max0(j-m,1)
 75 
            jk = max0(j-m,1)
77 
	    mu = max0(m+2-j,1)
 76 
            mu = max0(m+2-j,1)
78 
	    if (m .lt. mu) go to 20
 77 
            if (m .lt. mu) go to 20
79 
	    do 10 k = mu, m
 78 
            do 10 k = mu, m
80 
 
 79 
 
81 
	       t = abd(k,j) - ddot(k-mu,abd(ik,jk),1,abd(mu,j),1)
 80 
               t = abd(k,j) - ddot(k-mu,abd(ik,jk),1,abd(mu,j),1)
82 
	       t = t/abd(m+1,jk)
 81 
               t = t/abd(m+1,jk)
83 
	       abd(k,j) = t
 82 
               abd(k,j) = t
84 
	       s = s + t*t
 83 
               s = s + t*t
85 
	       ik = ik - 1
 84 
               ik = ik - 1
86 
	       jk = jk + 1
 85 
               jk = jk + 1
87 
   10	    continue
 86 
   10       continue
88 
   20	    continue
 87 
   20       continue
89 
 
 88 
 
90 
	    s = abd(m+1,j) - s
 89 
            s = abd(m+1,j) - s
91 
c     ......exit
 90 
c     ......exit
92 
	    if (s .le. 0.0e0) go to 40
 91 
            if (s .le. 0.0e0) go to 40
93 
 
 92 
 
94 
	    abd(m+1,j) = sqrt(s)
 93 
            abd(m+1,j) = sqrt(s)
95 
 
 94 
 
96 
   30	 continue
 95 
   30    continue
97 
	 info = 0
 96 
         info = 0
98 
   40 continue
 97 
   40 continue
99 
      return
 98 
      return
100 
      end
 99 
      end
101 
 
 100