1 /*
2 * Copyright (c) 2011, Michael Lehn
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1) Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2) Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided with the
15 * distribution.
16 * 3) Neither the name of the FLENS development group nor the names of
17 * its contributors may be used to endorse or promote products derived
18 * from this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 /* Based on
34 *
35 SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO )
36 *
37 * -- LAPACK routine (version 3.2) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * November 2006
41 */
42
43 #ifndef FLENS_LAPACK_GESV_TF2_TCC
44 #define FLENS_LAPACK_GESV_TF2_TCC 1
45
46 #include <algorithm>
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53
54 template <typename MA, typename VP>
55 typename GeMatrix<MA>::IndexType
56 tf2_generic(GeMatrix<MA> &A, DenseVector<VP> &piv)
57 {
58 using lapack::lamch;
59 using std::abs;
60 using std::min;
61
62 typedef typename GeMatrix<MA>::IndexType IndexType;
63 typedef typename GeMatrix<MA>::ElementType T;
64
65 const Underscore<IndexType> _;
66
67 const IndexType m = A.numRows();
68 const IndexType n = A.numCols();
69
70 IndexType info = 0;
71
72 //
73 // Quick return if possible
74 //
75 if ((m==0) || (n==0)) {
76 return info;
77 }
78 //
79 // Compute machine safe minimum
80 //
81 T safeMin = lamch<T>(SafeMin);
82
83 for (IndexType j=1; j<=min(m,n); ++j) {
84 //
85 // Row range of current submatrix A(j:M, j:N)
86 //
87 const auto rows = _(j, m);
88 //
89 // Row and column range of trailing submatrix A(j+1:M, j+1:N)
90 //
91 const auto _rows = _(j+1, m);
92 const auto _cols = _(j+1, n);
93 //
94 // Find pivot and test for singularity.
95 //
96 IndexType jp = j - 1 + blas::iamax(A(rows,j));
97 piv(j) = jp;
98 if (A(jp, j)!=T(0)) {
99 //
100 // Apply the interchange to columns 1:N.
101 //
102 if (j!=jp) {
103 blas::swap(A(j,_), A(jp,_));
104 }
105 //
106 // Compute elements J+1:M of J-th column.
107 //
108 if (j<m) {
109 if (abs(A(j,j))>=safeMin) {
110 blas::scal(T(1)/A(j, j), A(_rows,j));
111 } else {
112 for (IndexType i=1; i<=m-j; ++i) {
113 A(j+i,j) = A(j+i,j)/A(j,j);
114 }
115 }
116 }
117 } else {
118 if (info==0) {
119 info = j;
120 }
121 }
122 //
123 // Update trailing submatrix A(j+1:M, j+1:N)
124 //
125 blas::r(T(-1), A(_rows,j), A(j,_cols), A(_rows,_cols));
126 }
127
128 return info;
129 }
130
131 //== interface for native lapack ===============================================
132
133 #ifdef CHECK_CXXLAPACK
134
135 template <typename MA, typename VP>
136 typename GeMatrix<MA>::IndexType
137 tf2_native(GeMatrix<MA> &_A, DenseVector<VP> &piv)
138 {
139 typedef typename GeMatrix<MA>::ElementType T;
140
141 const INTEGER M = _A.numRows();
142 const INTEGER N = _A.numCols();
143 T *A = _A.data();
144 const INTEGER LDA = _A.leadingDimension();
145 INTEGER *IPIV = piv.data();
146 INTEGER INFO;
147
148 if (IsSame<T, DOUBLE>::value) {
149 LAPACK_IMPL(dgetf2)(&M, &N, A, &LDA, IPIV, &INFO);
150 } else {
151 ASSERT(0);
152 }
153 return INFO;
154 }
155
156 #endif // CHECK_CXXLAPACK
157
158 //== public interface ==========================================================
159
160 template <typename MA, typename VP>
161 typename GeMatrix<MA>::IndexType
162 tf2(GeMatrix<MA> &A, DenseVector<VP> &piv)
163 {
164 typedef typename GeMatrix<MA>::IndexType IndexType;
165
166 //
167 // Test the input parameters
168 //
169 ASSERT(A.firstRow()==1);
170 ASSERT(A.firstCol()==1);
171 ASSERT((piv.inc()>0 && piv.firstIndex()==1)
172 || (piv.inc()<0 && piv.firstIndex()==A.numRows()));
173
174 # ifdef CHECK_CXXLAPACK
175 //
176 // Make copies of output arguments
177 //
178 typename GeMatrix<MA>::NoView _A = A;
179 typename DenseVector<VP>::NoView _piv = piv;
180 # endif
181
182 //
183 // Call implementation
184 //
185 IndexType info = tf2_generic(A, piv);
186
187 # ifdef CHECK_CXXLAPACK
188 //
189 // Compare results
190 //
191 IndexType _info = tf2_native(_A, _piv);
192
193 bool failed = false;
194 if (! isIdentical(A, _A, " A", "_A")) {
195 std::cerr << "CXXLAPACK: A = " << A << std::endl;
196 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
197 failed = true;
198 }
199
200 if (! isIdentical(piv, _piv, " piv", "_piv")) {
201 std::cerr << "CXXLAPACK: piv = " << piv << std::endl;
202 std::cerr << "F77LAPACK: _piv = " << _piv << std::endl;
203 failed = true;
204 }
205
206 if (! isIdentical(info, _info, " info", "_info")) {
207 std::cerr << "CXXLAPACK: info = " << info << std::endl;
208 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
209 failed = true;
210 }
211
212 if (failed) {
213 ASSERT(0);
214 }
215 # endif
216
217 return info;
218 }
219
220 //-- forwarding ----------------------------------------------------------------
221 template <typename MA, typename VP>
222 typename MA::IndexType
223 tf2(MA &&A, VP &&piv)
224 {
225 typedef typename MA::IndexType IndexType;
226
227 CHECKPOINT_ENTER;
228 IndexType info = tf2(A, piv);
229 CHECKPOINT_LEAVE;
230
231 return info;
232 }
233
234 } } // namespace lapack, flens
235
236 #endif // FLENS_LAPACK_GESV_TF2_TCC
2 * Copyright (c) 2011, Michael Lehn
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1) Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2) Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided with the
15 * distribution.
16 * 3) Neither the name of the FLENS development group nor the names of
17 * its contributors may be used to endorse or promote products derived
18 * from this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 /* Based on
34 *
35 SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO )
36 *
37 * -- LAPACK routine (version 3.2) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * November 2006
41 */
42
43 #ifndef FLENS_LAPACK_GESV_TF2_TCC
44 #define FLENS_LAPACK_GESV_TF2_TCC 1
45
46 #include <algorithm>
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53
54 template <typename MA, typename VP>
55 typename GeMatrix<MA>::IndexType
56 tf2_generic(GeMatrix<MA> &A, DenseVector<VP> &piv)
57 {
58 using lapack::lamch;
59 using std::abs;
60 using std::min;
61
62 typedef typename GeMatrix<MA>::IndexType IndexType;
63 typedef typename GeMatrix<MA>::ElementType T;
64
65 const Underscore<IndexType> _;
66
67 const IndexType m = A.numRows();
68 const IndexType n = A.numCols();
69
70 IndexType info = 0;
71
72 //
73 // Quick return if possible
74 //
75 if ((m==0) || (n==0)) {
76 return info;
77 }
78 //
79 // Compute machine safe minimum
80 //
81 T safeMin = lamch<T>(SafeMin);
82
83 for (IndexType j=1; j<=min(m,n); ++j) {
84 //
85 // Row range of current submatrix A(j:M, j:N)
86 //
87 const auto rows = _(j, m);
88 //
89 // Row and column range of trailing submatrix A(j+1:M, j+1:N)
90 //
91 const auto _rows = _(j+1, m);
92 const auto _cols = _(j+1, n);
93 //
94 // Find pivot and test for singularity.
95 //
96 IndexType jp = j - 1 + blas::iamax(A(rows,j));
97 piv(j) = jp;
98 if (A(jp, j)!=T(0)) {
99 //
100 // Apply the interchange to columns 1:N.
101 //
102 if (j!=jp) {
103 blas::swap(A(j,_), A(jp,_));
104 }
105 //
106 // Compute elements J+1:M of J-th column.
107 //
108 if (j<m) {
109 if (abs(A(j,j))>=safeMin) {
110 blas::scal(T(1)/A(j, j), A(_rows,j));
111 } else {
112 for (IndexType i=1; i<=m-j; ++i) {
113 A(j+i,j) = A(j+i,j)/A(j,j);
114 }
115 }
116 }
117 } else {
118 if (info==0) {
119 info = j;
120 }
121 }
122 //
123 // Update trailing submatrix A(j+1:M, j+1:N)
124 //
125 blas::r(T(-1), A(_rows,j), A(j,_cols), A(_rows,_cols));
126 }
127
128 return info;
129 }
130
131 //== interface for native lapack ===============================================
132
133 #ifdef CHECK_CXXLAPACK
134
135 template <typename MA, typename VP>
136 typename GeMatrix<MA>::IndexType
137 tf2_native(GeMatrix<MA> &_A, DenseVector<VP> &piv)
138 {
139 typedef typename GeMatrix<MA>::ElementType T;
140
141 const INTEGER M = _A.numRows();
142 const INTEGER N = _A.numCols();
143 T *A = _A.data();
144 const INTEGER LDA = _A.leadingDimension();
145 INTEGER *IPIV = piv.data();
146 INTEGER INFO;
147
148 if (IsSame<T, DOUBLE>::value) {
149 LAPACK_IMPL(dgetf2)(&M, &N, A, &LDA, IPIV, &INFO);
150 } else {
151 ASSERT(0);
152 }
153 return INFO;
154 }
155
156 #endif // CHECK_CXXLAPACK
157
158 //== public interface ==========================================================
159
160 template <typename MA, typename VP>
161 typename GeMatrix<MA>::IndexType
162 tf2(GeMatrix<MA> &A, DenseVector<VP> &piv)
163 {
164 typedef typename GeMatrix<MA>::IndexType IndexType;
165
166 //
167 // Test the input parameters
168 //
169 ASSERT(A.firstRow()==1);
170 ASSERT(A.firstCol()==1);
171 ASSERT((piv.inc()>0 && piv.firstIndex()==1)
172 || (piv.inc()<0 && piv.firstIndex()==A.numRows()));
173
174 # ifdef CHECK_CXXLAPACK
175 //
176 // Make copies of output arguments
177 //
178 typename GeMatrix<MA>::NoView _A = A;
179 typename DenseVector<VP>::NoView _piv = piv;
180 # endif
181
182 //
183 // Call implementation
184 //
185 IndexType info = tf2_generic(A, piv);
186
187 # ifdef CHECK_CXXLAPACK
188 //
189 // Compare results
190 //
191 IndexType _info = tf2_native(_A, _piv);
192
193 bool failed = false;
194 if (! isIdentical(A, _A, " A", "_A")) {
195 std::cerr << "CXXLAPACK: A = " << A << std::endl;
196 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
197 failed = true;
198 }
199
200 if (! isIdentical(piv, _piv, " piv", "_piv")) {
201 std::cerr << "CXXLAPACK: piv = " << piv << std::endl;
202 std::cerr << "F77LAPACK: _piv = " << _piv << std::endl;
203 failed = true;
204 }
205
206 if (! isIdentical(info, _info, " info", "_info")) {
207 std::cerr << "CXXLAPACK: info = " << info << std::endl;
208 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
209 failed = true;
210 }
211
212 if (failed) {
213 ASSERT(0);
214 }
215 # endif
216
217 return info;
218 }
219
220 //-- forwarding ----------------------------------------------------------------
221 template <typename MA, typename VP>
222 typename MA::IndexType
223 tf2(MA &&A, VP &&piv)
224 {
225 typedef typename MA::IndexType IndexType;
226
227 CHECKPOINT_ENTER;
228 IndexType info = tf2(A, piv);
229 CHECKPOINT_LEAVE;
230
231 return info;
232 }
233
234 } } // namespace lapack, flens
235
236 #endif // FLENS_LAPACK_GESV_TF2_TCC