scid.ports.minpack.qform 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.qform;
9 
10 
11 private import std.algorithm: min;
12 
13 
14 
15 
16 /** This subroutine proceeds from the computed QR factorization of
17     an m by n matrix a to accumulate the m by m orthogonal matrix
18     Q from its factored form.
19 
20     Params:
21         m = a positive integer input variable set to the number
22             of rows of a and the order of q.
23         n = a positive integer input variable set to the number
24             of columns of a.
25         q = an array of length m^2. On input the full lower trapezoid in
26             the first min(m,n) columns of Q contains the factored form.
27             on output Q has been accumulated into a square matrix.
28         ldq = a positive integer input variable not less than m
29             which specifies the leading dimension of the array q.
30         wa = a work array of length m.
31 */
32 void qform(Real)(size_t m, size_t n, Real* q, size_t ldq, Real* wa)
33 {
34     size_t i, j, k, l;
35 
36     enum : Real
37     {
38         one = 1.0,
39         zero = 0.0
40     }
41 
42 
43     // Zero out upper triangle of Q in the first min(m,n) columns.
44     size_t minmn = min(m, n);
45     if (minmn >= 2)
46     {
47         for (j=1; j<minmn; j++)
48         {
49             for (i=0; i<j; i++)
50                 q[i+j*ldq] = zero;
51         }
52     }
53 
54 
55     // Initialize remaining columns to those of the identity matrix.
56     if (m > n)
57     {
58         for (j=n; j<m; j++)
59         {
60             for (i=0; i<m; i++)
61                 q[i+j*ldq] = zero;
62             q[j+j*ldq] = one;
63         }
64     }
65 
66 
67     // Accumulate Q from its factored form.
68     Real sum, temp;
69     for (l=0; l<minmn; l++)
70     {
71         k = minmn - l - 1;
72         for (i=k; i<m; i++)
73         {
74             wa[i] = q[i+k*ldq];
75             q[i+k*ldq] = zero;
76         }
77         q[k+k*ldq] = one;
78 
79         if (wa[k] != zero)
80         {
81             for (j=k; j<m; j++)
82             {
83                 sum = zero;
84                 for (i=k; i<m; i++)
85                     sum += q[i+j*ldq]*wa[i];
86                 temp = sum/wa[k];
87                 for (i=k; i<m; i++)
88                     q[i+j*ldq] = q[i+j*ldq] - temp*wa[i];
89             }
90         }
91     }
92 }
93 
94 
95 unittest { alias qform!(real) rqform; }