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 T Zero(0), One(1);
68
69 const IndexType m = A.numRows();
70 const IndexType n = A.numCols();
71
72 IndexType info = 0;
73
74 //
75 // Quick return if possible
76 //
77 if (m==0 || n==0) {
78 return info;
79 }
80 //
81 // Compute machine safe minimum
82 //
83 auto safeMin = lamch<T>(SafeMin);
84
85 for (IndexType j=1; j<=min(m,n); ++j) {
86 //
87 // Row range of current submatrix A(j:M, j:N)
88 //
89 const auto rows = _(j, m);
90 //
91 // Row and column range of trailing submatrix A(j+1:M, j+1:N)
92 //
93 const auto _rows = _(j+1, m);
94 const auto _cols = _(j+1, n);
95 //
96 // Find pivot and test for singularity.
97 //
98 IndexType jp = j - 1 + blas::iamax(A(rows,j));
99 piv(j) = jp;
100 if (A(jp,j)!=Zero) {
101 //
102 // Apply the interchange to columns 1:N.
103 //
104 if (j!=jp) {
105 blas::swap(A(j,_), A(jp,_));
106 }
107 //
108 // Compute elements J+1:M of J-th column.
109 //
110 if (j<m) {
111 if (abs(A(j,j))>=safeMin) {
112 blas::scal(One/A(j, j), A(_rows,j));
113 } else {
114 for (IndexType i=1; i<=m-j; ++i) {
115 A(j+i,j) /= A(j,j);
116 }
117 }
118 }
119 } else {
120 if (info==0) {
121 info = j;
122 }
123 }
124 //
125 // Update trailing submatrix A(j+1:M, j+1:N)
126 //
127 if (j<min(m,n)) {
128 blas::ru(-One, A(_rows,j), A(j,_cols), A(_rows,_cols));
129 }
130 }
131
132 return info;
133 }
134
135 //== interface for native lapack ===============================================
136
137 #ifdef CHECK_CXXLAPACK
138
139 template <typename MA, typename VP>
140 typename GeMatrix<MA>::IndexType
141 tf2_native(GeMatrix<MA> &A, DenseVector<VP> &piv)
142 {
143 return cxxlapack::getf2<INTEGER>(A.numRows(), A.numCols(),
144 A.data(), A.leadingDimension(),
145 piv.data());
146 }
147
148 #endif // CHECK_CXXLAPACK
149
150 //== public interface ==========================================================
151
152 template <typename MA, typename VP>
153 typename GeMatrix<MA>::IndexType
154 tf2(GeMatrix<MA> &A, DenseVector<VP> &piv)
155 {
156 typedef typename GeMatrix<MA>::IndexType IndexType;
157
158 //
159 // Test the input parameters
160 //
161 ASSERT(A.firstRow()==1);
162 ASSERT(A.firstCol()==1);
163 ASSERT((piv.inc()>0 && piv.firstIndex()==1)
164 || (piv.inc()<0 && piv.firstIndex()==A.numRows()));
165
166 # ifdef CHECK_CXXLAPACK
167 //
168 // Make copies of output arguments
169 //
170 typename GeMatrix<MA>::NoView A_org = A;
171 typename DenseVector<VP>::NoView piv_org = piv;
172 # endif
173
174 //
175 // Call implementation
176 //
177 std::cerr << "enter" << std::endl;
178 IndexType info = tf2_generic(A, piv);
179 std::cerr << "leave" << std::endl;
180
181 # ifdef CHECK_CXXLAPACK
182 //
183 // Make copies of results
184 //
185 typename GeMatrix<MA>::NoView A_generic = A;
186 typename DenseVector<VP>::NoView piv_generic = piv;
187 //
188 // Restore output arguments
189 //
190 A = A_org;
191 piv = piv_org;
192 //
193 // Compare results
194 //
195 IndexType _info = tf2_native(A, piv);
196
197 bool failed = false;
198 if (! isIdentical(A_generic, A, "A_generic", "A")) {
199 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
200 std::cerr << "F77LAPACK: A = " << A << std::endl;
201 failed = true;
202 }
203
204 if (! isIdentical(piv_generic, piv, "piv_generic", "piv")) {
205 std::cerr << "CXXLAPACK: piv_generic = " << piv_generic << std::endl;
206 std::cerr << "F77LAPACK: piv = " << piv << std::endl;
207 failed = true;
208 }
209
210 if (! isIdentical(info, _info, " info", "_info")) {
211 std::cerr << "CXXLAPACK: info = " << info << std::endl;
212 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
213 failed = true;
214 }
215
216 if (failed) {
217 std::cerr << "A_org = " << A_org << std::endl;
218 ASSERT(0);
219 }
220 # endif
221
222 return info;
223 }
224
225 //-- forwarding ----------------------------------------------------------------
226 template <typename MA, typename VP>
227 typename MA::IndexType
228 tf2(MA &&A, VP &&piv)
229 {
230 typedef typename MA::IndexType IndexType;
231
232 CHECKPOINT_ENTER;
233 IndexType info = tf2(A, piv);
234 CHECKPOINT_LEAVE;
235
236 return info;
237 }
238
239 } } // namespace lapack, flens
240
241 #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 T Zero(0), One(1);
68
69 const IndexType m = A.numRows();
70 const IndexType n = A.numCols();
71
72 IndexType info = 0;
73
74 //
75 // Quick return if possible
76 //
77 if (m==0 || n==0) {
78 return info;
79 }
80 //
81 // Compute machine safe minimum
82 //
83 auto safeMin = lamch<T>(SafeMin);
84
85 for (IndexType j=1; j<=min(m,n); ++j) {
86 //
87 // Row range of current submatrix A(j:M, j:N)
88 //
89 const auto rows = _(j, m);
90 //
91 // Row and column range of trailing submatrix A(j+1:M, j+1:N)
92 //
93 const auto _rows = _(j+1, m);
94 const auto _cols = _(j+1, n);
95 //
96 // Find pivot and test for singularity.
97 //
98 IndexType jp = j - 1 + blas::iamax(A(rows,j));
99 piv(j) = jp;
100 if (A(jp,j)!=Zero) {
101 //
102 // Apply the interchange to columns 1:N.
103 //
104 if (j!=jp) {
105 blas::swap(A(j,_), A(jp,_));
106 }
107 //
108 // Compute elements J+1:M of J-th column.
109 //
110 if (j<m) {
111 if (abs(A(j,j))>=safeMin) {
112 blas::scal(One/A(j, j), A(_rows,j));
113 } else {
114 for (IndexType i=1; i<=m-j; ++i) {
115 A(j+i,j) /= A(j,j);
116 }
117 }
118 }
119 } else {
120 if (info==0) {
121 info = j;
122 }
123 }
124 //
125 // Update trailing submatrix A(j+1:M, j+1:N)
126 //
127 if (j<min(m,n)) {
128 blas::ru(-One, A(_rows,j), A(j,_cols), A(_rows,_cols));
129 }
130 }
131
132 return info;
133 }
134
135 //== interface for native lapack ===============================================
136
137 #ifdef CHECK_CXXLAPACK
138
139 template <typename MA, typename VP>
140 typename GeMatrix<MA>::IndexType
141 tf2_native(GeMatrix<MA> &A, DenseVector<VP> &piv)
142 {
143 return cxxlapack::getf2<INTEGER>(A.numRows(), A.numCols(),
144 A.data(), A.leadingDimension(),
145 piv.data());
146 }
147
148 #endif // CHECK_CXXLAPACK
149
150 //== public interface ==========================================================
151
152 template <typename MA, typename VP>
153 typename GeMatrix<MA>::IndexType
154 tf2(GeMatrix<MA> &A, DenseVector<VP> &piv)
155 {
156 typedef typename GeMatrix<MA>::IndexType IndexType;
157
158 //
159 // Test the input parameters
160 //
161 ASSERT(A.firstRow()==1);
162 ASSERT(A.firstCol()==1);
163 ASSERT((piv.inc()>0 && piv.firstIndex()==1)
164 || (piv.inc()<0 && piv.firstIndex()==A.numRows()));
165
166 # ifdef CHECK_CXXLAPACK
167 //
168 // Make copies of output arguments
169 //
170 typename GeMatrix<MA>::NoView A_org = A;
171 typename DenseVector<VP>::NoView piv_org = piv;
172 # endif
173
174 //
175 // Call implementation
176 //
177 std::cerr << "enter" << std::endl;
178 IndexType info = tf2_generic(A, piv);
179 std::cerr << "leave" << std::endl;
180
181 # ifdef CHECK_CXXLAPACK
182 //
183 // Make copies of results
184 //
185 typename GeMatrix<MA>::NoView A_generic = A;
186 typename DenseVector<VP>::NoView piv_generic = piv;
187 //
188 // Restore output arguments
189 //
190 A = A_org;
191 piv = piv_org;
192 //
193 // Compare results
194 //
195 IndexType _info = tf2_native(A, piv);
196
197 bool failed = false;
198 if (! isIdentical(A_generic, A, "A_generic", "A")) {
199 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
200 std::cerr << "F77LAPACK: A = " << A << std::endl;
201 failed = true;
202 }
203
204 if (! isIdentical(piv_generic, piv, "piv_generic", "piv")) {
205 std::cerr << "CXXLAPACK: piv_generic = " << piv_generic << std::endl;
206 std::cerr << "F77LAPACK: piv = " << piv << std::endl;
207 failed = true;
208 }
209
210 if (! isIdentical(info, _info, " info", "_info")) {
211 std::cerr << "CXXLAPACK: info = " << info << std::endl;
212 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
213 failed = true;
214 }
215
216 if (failed) {
217 std::cerr << "A_org = " << A_org << std::endl;
218 ASSERT(0);
219 }
220 # endif
221
222 return info;
223 }
224
225 //-- forwarding ----------------------------------------------------------------
226 template <typename MA, typename VP>
227 typename MA::IndexType
228 tf2(MA &&A, VP &&piv)
229 {
230 typedef typename MA::IndexType IndexType;
231
232 CHECKPOINT_ENTER;
233 IndexType info = tf2(A, piv);
234 CHECKPOINT_LEAVE;
235
236 return info;
237 }
238
239 } } // namespace lapack, flens
240
241 #endif // FLENS_LAPACK_GESV_TF2_TCC