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 DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK,
36 $ INFO )
37 *
38 * -- LAPACK routine (version 3.3.1) --
39 * -- LAPACK is a software package provided by Univ. of Tennessee, --
40 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
41 * -- April 2011 --
42 */
43
44 #ifndef FLENS_LAPACK_AUX_CON_TCC
45 #define FLENS_LAPACK_AUX_CON_TCC 1
46
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53 template <typename MA, typename NORMA, typename RCOND,
54 typename VWORK, typename VIWORK>
55 void
56 con_generic(Norm norm,
57 const GeMatrix<MA> &A,
58 const NORMA &normA,
59 RCOND &rCond,
60 DenseVector<VWORK> &work,
61 DenseVector<VIWORK> &iwork)
62 {
63 using std::abs;
64
65 typedef typename GeMatrix<MA>::ElementType ElementType;
66 typedef typename GeMatrix<MA>::IndexType IndexType;
67
68 const ElementType Zero(0), One(1);
69
70 const Underscore<IndexType> _;
71
72 const IndexType n = A.numRows();
73
74 //
75 // Local Arrays
76 //
77 IndexType iSaveData[3] = {0, 0, 0};
78 DenseVectorView<IndexType>
79 iSave = typename DenseVectorView<IndexType>::Engine(3, iSaveData, 1);
80 //
81 // Quick return if possible
82 //
83 rCond = Zero;
84 if (n==0) {
85 rCond = One;
86 return;
87 } else if (normA==Zero) {
88 return;
89 }
90
91 const ElementType smallNum = lamch<ElementType>(SafeMin);
92 //
93 // Estimate the norm of inv(A).
94 //
95 ElementType normInvA = Zero;
96 bool normIn = false;
97
98 IndexType kase, kase1;
99
100 if (norm==OneNorm) {
101 kase1 = 1;
102 } else {
103 kase1 = 2;
104 }
105 kase = 0;
106
107 auto work1 = work(_(1,n));
108 auto work2 = work(_(n+1,2*n));
109 auto work3 = work(_(2*n+1,3*n));
110 auto work4 = work(_(3*n+1,4*n));
111
112 bool computeRCond = true;
113
114 while (true) {
115 lacn2(work2, work1, iwork, normInvA, kase, iSave);
116 if (kase==0) {
117 break;
118 }
119
120 ElementType sl, su;
121 if (kase==kase1) {
122 //
123 // Multiply by inv(L).
124 //
125 latrs(NoTrans, normIn, A.lowerUnit(), work1, sl, work3);
126 //
127 // Multiply by inv(U).
128 //
129 latrs(NoTrans, normIn, A.upper(), work1, su, work4);
130 } else {
131 //
132 // Multiply by inv(U**T).
133 //
134 latrs(Trans, normIn, A.upper(), work1, su, work4);
135 //
136 // Multiply by inv(L**T).
137 //
138 latrs(Trans, normIn, A.lowerUnit(), work1, sl, work3);
139 }
140 //
141 // Divide X by 1/(SL*SU) if doing so will not cause overflow.
142 //
143 const ElementType scale = sl*su;
144 normIn = true;
145 if (scale!=One) {
146 const IndexType ix = blas::iamax(work1);
147 if (scale<abs(work1(ix))*smallNum || scale==Zero) {
148 computeRCond = false;
149 break;
150 }
151 rscl(scale, work1);
152 }
153 }
154 //
155 // Compute the estimate of the reciprocal condition number.
156 //
157 if (computeRCond) {
158 if (normInvA!=Zero) {
159 rCond = (One/normInvA) / normA;
160 }
161 }
162 }
163
164 //== interface for native lapack ===============================================
165
166 #ifdef CHECK_CXXLAPACK
167
168 template <typename MA, typename NORMA, typename RCOND,
169 typename VWORK, typename VIWORK>
170 void
171 con_native(Norm norm,
172 const GeMatrix<MA> &A,
173 const NORMA &normA,
174 RCOND &rCond,
175 DenseVector<VWORK> &work,
176 DenseVector<VIWORK> &iwork)
177 {
178 typedef typename GeMatrix<MA>::ElementType T;
179
180 const char NORM = char(norm);
181 const INTEGER N = A.numRows();
182 const INTEGER LDA = A.leadingDimension();
183 const T ANORM = normA;
184 INTEGER INFO;
185
186 if (IsSame<T,double>::value) {
187 LAPACK_IMPL(dgecon)(&NORM,
188 &N,
189 A.data(),
190 &LDA,
191 &ANORM,
192 &rCond,
193 work.data(),
194 iwork.data(),
195 &INFO);
196 } else {
197 ASSERT(0);
198 }
199 ASSERT(INFO>=0);
200 }
201
202 #endif // CHECK_CXXLAPACK
203
204 //== public interface ==========================================================
205 template <typename MA, typename NORMA, typename RCOND,
206 typename VWORK, typename VIWORK>
207 void
208 con(Norm norm,
209 const GeMatrix<MA> &A,
210 const NORMA &normA,
211 RCOND &rCond,
212 DenseVector<VWORK> &work,
213 DenseVector<VIWORK> &iwork)
214 {
215 typedef typename GeMatrix<MA>::IndexType IndexType;
216 //
217 // Test the input parameters
218 //
219 # ifndef NDEBUG
220 ASSERT(norm==InfinityNorm || norm==OneNorm);
221 ASSERT(A.firstRow()==1);
222 ASSERT(A.firstCol()==1);
223 ASSERT(A.numRows()==A.numCols());
224
225 const IndexType n = A.numRows();
226
227 ASSERT(work.firstIndex()==1);
228 ASSERT(work.length()==4*n);
229
230 ASSERT(iwork.firstIndex()==1);
231 ASSERT(iwork.length()==n);
232 # endif
233
234 # ifdef CHECK_CXXLAPACK
235 //
236 // Make copies of output arguments
237 //
238 RCOND rCond_org = rCond;
239 typename DenseVector<VWORK>::NoView work_org = work;
240 typename DenseVector<VIWORK>::NoView iwork_org = iwork;
241 # endif
242
243 //
244 // Call implementation
245 //
246 con_generic(norm, A, normA, rCond, work, iwork);
247
248 # ifdef CHECK_CXXLAPACK
249 //
250 // Compare results
251 //
252 RCOND rCond_generic = rCond;
253 typename DenseVector<VWORK>::NoView work_generic = work;
254 typename DenseVector<VIWORK>::NoView iwork_generic = iwork;
255
256 rCond = rCond_org;
257 work = work_org;
258 iwork = iwork_org;
259
260 con_native(norm, A, normA, rCond, work, iwork);
261
262 bool failed = false;
263 if (! isIdentical(rCond_generic, rCond, "rCond_generic", "rCond")) {
264 std::cerr << "CXXLAPACK: rCond_generic = "
265 << rCond_generic << std::endl;
266 std::cerr << "F77LAPACK: rCond = " << rCond << std::endl;
267 failed = true;
268 }
269
270 if (! isIdentical(work_generic, work, "work_generic", "work")) {
271 std::cerr << "CXXLAPACK: work_generic = "
272 << work_generic << std::endl;
273 std::cerr << "F77LAPACK: work = " << work << std::endl;
274 failed = true;
275 }
276
277 if (! isIdentical(iwork_generic, iwork, "iwork_generic", "iwork")) {
278 std::cerr << "CXXLAPACK: iwork_generic = "
279 << iwork_generic << std::endl;
280 std::cerr << "F77LAPACK: iwork = " << iwork << std::endl;
281 failed = true;
282 }
283
284 if (failed) {
285 ASSERT(0);
286 } else {
287 // std::cerr << "passed: con" << std::endl;
288 }
289 # endif
290 }
291
292 //-- forwarding ----------------------------------------------------------------
293 template <typename MA, typename NORMA, typename RCOND,
294 typename VWORK, typename VIWORK>
295 void
296 con(Norm norm,
297 const MA &A,
298 const NORMA &normA,
299 RCOND &&rCond,
300 VWORK &&work,
301 VIWORK &&iwork)
302 {
303 CHECKPOINT_ENTER;
304 con(norm, A, normA, rCond, work, iwork);
305 CHECKPOINT_LEAVE;
306 }
307
308 } } // namespace lapack, flens
309
310 #endif // FLENS_LAPACK_AUX_CON_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 DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK,
36 $ INFO )
37 *
38 * -- LAPACK routine (version 3.3.1) --
39 * -- LAPACK is a software package provided by Univ. of Tennessee, --
40 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
41 * -- April 2011 --
42 */
43
44 #ifndef FLENS_LAPACK_AUX_CON_TCC
45 #define FLENS_LAPACK_AUX_CON_TCC 1
46
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53 template <typename MA, typename NORMA, typename RCOND,
54 typename VWORK, typename VIWORK>
55 void
56 con_generic(Norm norm,
57 const GeMatrix<MA> &A,
58 const NORMA &normA,
59 RCOND &rCond,
60 DenseVector<VWORK> &work,
61 DenseVector<VIWORK> &iwork)
62 {
63 using std::abs;
64
65 typedef typename GeMatrix<MA>::ElementType ElementType;
66 typedef typename GeMatrix<MA>::IndexType IndexType;
67
68 const ElementType Zero(0), One(1);
69
70 const Underscore<IndexType> _;
71
72 const IndexType n = A.numRows();
73
74 //
75 // Local Arrays
76 //
77 IndexType iSaveData[3] = {0, 0, 0};
78 DenseVectorView<IndexType>
79 iSave = typename DenseVectorView<IndexType>::Engine(3, iSaveData, 1);
80 //
81 // Quick return if possible
82 //
83 rCond = Zero;
84 if (n==0) {
85 rCond = One;
86 return;
87 } else if (normA==Zero) {
88 return;
89 }
90
91 const ElementType smallNum = lamch<ElementType>(SafeMin);
92 //
93 // Estimate the norm of inv(A).
94 //
95 ElementType normInvA = Zero;
96 bool normIn = false;
97
98 IndexType kase, kase1;
99
100 if (norm==OneNorm) {
101 kase1 = 1;
102 } else {
103 kase1 = 2;
104 }
105 kase = 0;
106
107 auto work1 = work(_(1,n));
108 auto work2 = work(_(n+1,2*n));
109 auto work3 = work(_(2*n+1,3*n));
110 auto work4 = work(_(3*n+1,4*n));
111
112 bool computeRCond = true;
113
114 while (true) {
115 lacn2(work2, work1, iwork, normInvA, kase, iSave);
116 if (kase==0) {
117 break;
118 }
119
120 ElementType sl, su;
121 if (kase==kase1) {
122 //
123 // Multiply by inv(L).
124 //
125 latrs(NoTrans, normIn, A.lowerUnit(), work1, sl, work3);
126 //
127 // Multiply by inv(U).
128 //
129 latrs(NoTrans, normIn, A.upper(), work1, su, work4);
130 } else {
131 //
132 // Multiply by inv(U**T).
133 //
134 latrs(Trans, normIn, A.upper(), work1, su, work4);
135 //
136 // Multiply by inv(L**T).
137 //
138 latrs(Trans, normIn, A.lowerUnit(), work1, sl, work3);
139 }
140 //
141 // Divide X by 1/(SL*SU) if doing so will not cause overflow.
142 //
143 const ElementType scale = sl*su;
144 normIn = true;
145 if (scale!=One) {
146 const IndexType ix = blas::iamax(work1);
147 if (scale<abs(work1(ix))*smallNum || scale==Zero) {
148 computeRCond = false;
149 break;
150 }
151 rscl(scale, work1);
152 }
153 }
154 //
155 // Compute the estimate of the reciprocal condition number.
156 //
157 if (computeRCond) {
158 if (normInvA!=Zero) {
159 rCond = (One/normInvA) / normA;
160 }
161 }
162 }
163
164 //== interface for native lapack ===============================================
165
166 #ifdef CHECK_CXXLAPACK
167
168 template <typename MA, typename NORMA, typename RCOND,
169 typename VWORK, typename VIWORK>
170 void
171 con_native(Norm norm,
172 const GeMatrix<MA> &A,
173 const NORMA &normA,
174 RCOND &rCond,
175 DenseVector<VWORK> &work,
176 DenseVector<VIWORK> &iwork)
177 {
178 typedef typename GeMatrix<MA>::ElementType T;
179
180 const char NORM = char(norm);
181 const INTEGER N = A.numRows();
182 const INTEGER LDA = A.leadingDimension();
183 const T ANORM = normA;
184 INTEGER INFO;
185
186 if (IsSame<T,double>::value) {
187 LAPACK_IMPL(dgecon)(&NORM,
188 &N,
189 A.data(),
190 &LDA,
191 &ANORM,
192 &rCond,
193 work.data(),
194 iwork.data(),
195 &INFO);
196 } else {
197 ASSERT(0);
198 }
199 ASSERT(INFO>=0);
200 }
201
202 #endif // CHECK_CXXLAPACK
203
204 //== public interface ==========================================================
205 template <typename MA, typename NORMA, typename RCOND,
206 typename VWORK, typename VIWORK>
207 void
208 con(Norm norm,
209 const GeMatrix<MA> &A,
210 const NORMA &normA,
211 RCOND &rCond,
212 DenseVector<VWORK> &work,
213 DenseVector<VIWORK> &iwork)
214 {
215 typedef typename GeMatrix<MA>::IndexType IndexType;
216 //
217 // Test the input parameters
218 //
219 # ifndef NDEBUG
220 ASSERT(norm==InfinityNorm || norm==OneNorm);
221 ASSERT(A.firstRow()==1);
222 ASSERT(A.firstCol()==1);
223 ASSERT(A.numRows()==A.numCols());
224
225 const IndexType n = A.numRows();
226
227 ASSERT(work.firstIndex()==1);
228 ASSERT(work.length()==4*n);
229
230 ASSERT(iwork.firstIndex()==1);
231 ASSERT(iwork.length()==n);
232 # endif
233
234 # ifdef CHECK_CXXLAPACK
235 //
236 // Make copies of output arguments
237 //
238 RCOND rCond_org = rCond;
239 typename DenseVector<VWORK>::NoView work_org = work;
240 typename DenseVector<VIWORK>::NoView iwork_org = iwork;
241 # endif
242
243 //
244 // Call implementation
245 //
246 con_generic(norm, A, normA, rCond, work, iwork);
247
248 # ifdef CHECK_CXXLAPACK
249 //
250 // Compare results
251 //
252 RCOND rCond_generic = rCond;
253 typename DenseVector<VWORK>::NoView work_generic = work;
254 typename DenseVector<VIWORK>::NoView iwork_generic = iwork;
255
256 rCond = rCond_org;
257 work = work_org;
258 iwork = iwork_org;
259
260 con_native(norm, A, normA, rCond, work, iwork);
261
262 bool failed = false;
263 if (! isIdentical(rCond_generic, rCond, "rCond_generic", "rCond")) {
264 std::cerr << "CXXLAPACK: rCond_generic = "
265 << rCond_generic << std::endl;
266 std::cerr << "F77LAPACK: rCond = " << rCond << std::endl;
267 failed = true;
268 }
269
270 if (! isIdentical(work_generic, work, "work_generic", "work")) {
271 std::cerr << "CXXLAPACK: work_generic = "
272 << work_generic << std::endl;
273 std::cerr << "F77LAPACK: work = " << work << std::endl;
274 failed = true;
275 }
276
277 if (! isIdentical(iwork_generic, iwork, "iwork_generic", "iwork")) {
278 std::cerr << "CXXLAPACK: iwork_generic = "
279 << iwork_generic << std::endl;
280 std::cerr << "F77LAPACK: iwork = " << iwork << std::endl;
281 failed = true;
282 }
283
284 if (failed) {
285 ASSERT(0);
286 } else {
287 // std::cerr << "passed: con" << std::endl;
288 }
289 # endif
290 }
291
292 //-- forwarding ----------------------------------------------------------------
293 template <typename MA, typename NORMA, typename RCOND,
294 typename VWORK, typename VIWORK>
295 void
296 con(Norm norm,
297 const MA &A,
298 const NORMA &normA,
299 RCOND &&rCond,
300 VWORK &&work,
301 VIWORK &&iwork)
302 {
303 CHECKPOINT_ENTER;
304 con(norm, A, normA, rCond, work, iwork);
305 CHECKPOINT_LEAVE;
306 }
307
308 } } // namespace lapack, flens
309
310 #endif // FLENS_LAPACK_AUX_CON_TCC