← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/optimize/minpack/r1mpyq.f

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/optimize/minpack/r1mpyq.f
 
1
      subroutine r1mpyq(m,n,a,lda,v,w)
 
2
      integer m,n,lda
 
3
      double precision a(lda,n),v(n),w(n)
 
4
c     **********
 
5
c
 
6
c     subroutine r1mpyq
 
7
c
 
8
c     given an m by n matrix a, this subroutine computes a*q where
 
9
c     q is the product of 2*(n - 1) transformations
 
10
c
 
11
c           gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
 
12
c
 
13
c     and gv(i), gw(i) are givens rotations in the (i,n) plane which
 
14
c     eliminate elements in the i-th and n-th planes, respectively.
 
15
c     q itself is not given, rather the information to recover the
 
16
c     gv, gw rotations is supplied.
 
17
c
 
18
c     the subroutine statement is
 
19
c
 
20
c       subroutine r1mpyq(m,n,a,lda,v,w)
 
21
c
 
22
c     where
 
23
c
 
24
c       m is a positive integer input variable set to the number
 
25
c         of rows of a.
 
26
c
 
27
c       n is a positive integer input variable set to the number
 
28
c         of columns of a.
 
29
c
 
30
c       a is an m by n array. on input a must contain the matrix
 
31
c         to be postmultiplied by the orthogonal matrix q
 
32
c         described above. on output a*q has replaced a.
 
33
c
 
34
c       lda is a positive integer input variable not less than m
 
35
c         which specifies the leading dimension of the array a.
 
36
c
 
37
c       v is an input array of length n. v(i) must contain the
 
38
c         information necessary to recover the givens rotation gv(i)
 
39
c         described above.
 
40
c
 
41
c       w is an input array of length n. w(i) must contain the
 
42
c         information necessary to recover the givens rotation gw(i)
 
43
c         described above.
 
44
c
 
45
c     subroutines called
 
46
c
 
47
c       fortran-supplied ... dabs,dsqrt
 
48
c
 
49
c     argonne national laboratory. minpack project. march 1980.
 
50
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
 
51
c
 
52
c     **********
 
53
      integer i,j,nmj,nm1
 
54
      double precision cos,one,sin,temp
 
55
      data one /1.0d0/
 
56
c
 
57
c     apply the first set of givens rotations to a.
 
58
c
 
59
      nm1 = n - 1
 
60
      if (nm1 .lt. 1) go to 50
 
61
      do 20 nmj = 1, nm1
 
62
         j = n - nmj
 
63
         if (dabs(v(j)) .gt. one) cos = one/v(j)
 
64
         if (dabs(v(j)) .gt. one) sin = dsqrt(one-cos**2)
 
65
         if (dabs(v(j)) .le. one) sin = v(j)
 
66
         if (dabs(v(j)) .le. one) cos = dsqrt(one-sin**2)
 
67
         do 10 i = 1, m
 
68
            temp = cos*a(i,j) - sin*a(i,n)
 
69
            a(i,n) = sin*a(i,j) + cos*a(i,n)
 
70
            a(i,j) = temp
 
71
   10       continue
 
72
   20    continue
 
73
c
 
74
c     apply the second set of givens rotations to a.
 
75
c
 
76
      do 40 j = 1, nm1
 
77
         if (dabs(w(j)) .gt. one) cos = one/w(j)
 
78
         if (dabs(w(j)) .gt. one) sin = dsqrt(one-cos**2)
 
79
         if (dabs(w(j)) .le. one) sin = w(j)
 
80
         if (dabs(w(j)) .le. one) cos = dsqrt(one-sin**2)
 
81
         do 30 i = 1, m
 
82
            temp = cos*a(i,j) + sin*a(i,n)
 
83
            a(i,n) = -sin*a(i,j) + cos*a(i,n)
 
84
            a(i,j) = temp
 
85
   30       continue
 
86
   40    continue
 
87
   50 continue
 
88
      return
 
89
c
 
90
c     last card of subroutine r1mpyq.
 
91
c
 
92
      end