scid.ports.minpack.r1mpyq source code

1 // This code has been mechanically translated from the original FORTRAN
2 // code at http://netlib.org/minpack.
3 
4 /** Authors:    Lars Tandle Kyllingstad
5     Copyright:  Copyright (c) 2009, Lars T. Kyllingstad. All rights reserved.
6     License:    Boost Licence 1.0
7 */
8 module scid.ports.minpack.r1mpyq;
9 
10 
11 private import std.math: abs, sqrt;
12 
13 
14 
15 
16 /** Given an m by n matrix A, this subroutine computes AQ where
17     Q is the product of 2*(n-1) transformations
18     ---
19     gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
20     ---
21     and gv(i), gw(i) are Givens rotations in the (i,n) plane which
22     eliminate elements in the i-th and n-th planes, respectively.
23     Q itself is not given, rather the information to recover the
24     gv, gw rotations is supplied.
25 
26 
27     Params:
28         m = a positive integer input variable set to the number
29             of rows of A.
30         
31         n = a positive integer input variable set to the number
32             of columns of A.
33         
34         a = an array of length m^2. On input a must contain the matrix
35             A to be postmultiplied by the orthogonal matrix Q
36             described above. On output AQ has replaced A.
37             the first min(m,n) columns of Q contains the factored form.
38             on output Q has been accumulated into a square matrix.
39         
40         lda = a positive integer input variable not less than m
41             which specifies the leading dimension of the array a.
42         
43         v = an input array of length n. v(i) must contain the
44             information necessary to recover the Givens rotation gv(i)
45             described above.
46 
47         v = an input array of length n. v(i) must contain the
48             information necessary to recover the Givens rotation gw(i)
49             described above.
50 */
51 void r1mpyq(Real)(size_t m, size_t n, Real* a, size_t lda, Real* v, Real* w)
52 {
53     size_t i, j, ij, in_, nmj;
54     Real cos, sin, temp;
55     
56     enum  Real one = 1.0;
57 
58     
59     // Apply the first set of Givens rotations to A.
60     if (n > 1)
61     {
62         size_t nm1 = n - 1;
63         for (nmj=2; nmj<=n; nmj++)
64         {
65             j = n - nmj;
66             if (abs(v[j]) > one)
67             {
68                 cos = one/v[j];
69                 sin = sqrt(one-cos*cos);
70             }
71             else
72             {
73                 sin = v[j];
74                 cos = sqrt(one-sin*sin);
75             }
76 
77             for (i=0; i<m; i++)
78             {
79                 ij = i + j*lda;
80                 in_ = i + nm1*lda;
81                 temp = cos*a[ij] - sin*a[in_];
82                 a[in_] = sin*a[ij] + cos*a[in_];
83                 a[ij] = temp;
84             }
85         }
86 
87 
88         // Apply the second set of Givens rotations to A.
89         for (j=0; j<nm1; j++)
90         {
91             if (abs(w[j]) > one)
92             {
93                 cos = one/w[j];
94                 sin = sqrt(one-cos*cos);
95             }
96             else
97             {
98                 sin = w[j];
99                 cos = sqrt(one-sin*sin);
100             }
101 
102             for (i=0; i<m; i++)
103             {
104                 ij = i + j*lda;
105                 in_ = i + nm1*lda;
106                 temp = cos*a[ij] + sin*a[in_];
107                 a[in_] = -sin*a[ij] + cos*a[in_];
108                 a[ij] = temp;
109             }
110         }
111     } 
112 }
113 
114 
115 unittest { alias r1mpyq!(real) rr1mpyq; }