← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/fpsysy.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/interpolate/fitpack/fpsysy.f
1
 
      subroutine fpsysy(a,n,g)
2
 
c subroutine fpsysy solves a linear n x n symmetric system
3
 
c    (a) * (b) = (g)
4
 
c on input, vector g contains the right hand side ; on output it will
5
 
c contain the solution (b).
6
 
c  ..
7
 
c  ..scalar arguments..
8
 
      integer n
9
 
c  ..array arguments..
10
 
      real*8 a(6,6),g(6)
11
 
c  ..local scalars..
12
 
      real*8 fac
13
 
      integer i,i1,j,k
14
 
c  ..
15
 
      g(1) = g(1)/a(1,1)
16
 
      if(n.eq.1) return
17
 
c  decomposition of the symmetric matrix (a) = (l) * (d) *(l)'
18
 
c  with (l) a unit lower triangular matrix and (d) a diagonal
19
 
c  matrix
20
 
      do 10 k=2,n
21
 
         a(k,1) = a(k,1)/a(1,1)
22
 
  10  continue
23
 
      do 40 i=2,n
24
 
         i1 = i-1
25
 
         do 30 k=i,n
26
 
            fac = a(k,i)
27
 
            do 20 j=1,i1
28
 
               fac = fac-a(j,j)*a(k,j)*a(i,j)
29
 
  20        continue
30
 
            a(k,i) = fac
31
 
            if(k.gt.i) a(k,i) = fac/a(i,i)
32
 
  30     continue
33
 
  40  continue
34
 
c  solve the system (l)*(d)*(l)'*(b) = (g).
35
 
c  first step : solve (l)*(d)*(c) = (g).
36
 
      do 60 i=2,n
37
 
         i1 = i-1
38
 
         fac = g(i)
39
 
         do 50 j=1,i1
40
 
            fac = fac-g(j)*a(j,j)*a(i,j)
41
 
  50     continue
42
 
         g(i) = fac/a(i,i)
43
 
  60  continue
44
 
c  second step : solve (l)'*(b) = (c)
45
 
      i = n
46
 
      do 80 j=2,n
47
 
         i1 = i
48
 
         i = i-1
49
 
         fac = g(i)
50
 
         do 70 k=i1,n
51
 
            fac = fac-g(k)*a(k,i)
52
 
  70     continue
53
 
         g(i) = fac
54
 
  80  continue
55
 
      return
56
 
      end