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 DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
36 *
37 * -- LAPACK driver 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_SV_TCC
44 #define FLENS_LAPACK_GESV_SV_TCC 1
45
46 #include <flens/blas/blas.h>
47 #include <flens/lapack/lapack.h>
48
49 namespace flens { namespace lapack {
50
51 //== generic lapack implementation =============================================
52 template <typename MA, typename VP, typename MB>
53 typename GeMatrix<MA>::IndexType
54 sv_generic(GeMatrix<MA> &A, DenseVector<VP> &piv, GeMatrix<MB> &B)
55 {
56 typedef typename GeMatrix<MA>::IndexType IndexType;
57
58 IndexType info = 0;
59
60 //
61 // Compute the LU factorization of A.
62 //
63 info = trf(A, piv);
64 if (info==0) {
65 //
66 // Solve the system A*X = B, overwriting B with X.
67 //
68 trs(NoTrans, A, piv, B);
69 }
70 return info;
71 }
72
73 //== interface for native lapack ===============================================
74
75 #ifdef CHECK_CXXLAPACK
76 template <typename MA, typename VP, typename MB>
77 typename GeMatrix<MA>::IndexType
78 sv_native(GeMatrix<MA> &A, DenseVector<VP> &piv, GeMatrix<MB> &B)
79 {
80 typedef typename GeMatrix<MA>::ElementType ElementType;
81 const INTEGER N = A.numRows();
82 const INTEGER NRHS = B.numCols();
83 const INTEGER LDA = A.leadingDimension();
84 const INTEGER LDB = B.leadingDimension();
85 INTEGER INFO;
86
87 if (IsSame<ElementType, DOUBLE>::value) {
88 LAPACK_IMPL(dgesv)(&N,
89 &NRHS,
90 A.data(),
91 &LDA,
92 piv.data(),
93 B.data(),
94 &LDB,
95 &INFO);
96 } else {
97 ASSERT(0);
98 }
99 return INFO;
100 }
101
102 #endif // CHECK_CXXLAPACK
103
104 //== public interface ==========================================================
105
106 template <typename MA, typename VP, typename MB>
107 typename GeMatrix<MA>::IndexType
108 sv(GeMatrix<MA> &A, DenseVector<VP> &piv, GeMatrix<MB> &B)
109 {
110 typedef typename GeMatrix<MA>::IndexType IndexType;
111 //
112 // Test the input parameters
113 //
114 # ifndef NDEBUG
115 ASSERT(A.firstRow()==1);
116 ASSERT(A.firstCol()==1);
117 ASSERT(A.numRows()==A.numCols());
118 ASSERT((piv.inc()>0 && piv.firstIndex()==1)
119 || (piv.inc()<0 && piv.firstIndex()==A.numRows()));
120 ASSERT(B.firstRow()==1);
121 ASSERT(B.firstCol()==1);
122 ASSERT(B.numRows()==A.numRows());
123 # endif
124 //
125 // Make copies of output arguments
126 //
127 typename GeMatrix<MA>::NoView A_org = A;
128 typename DenseVector<VP>::NoView piv_org = piv;
129 typename GeMatrix<MB>::NoView B_org = B;
130 //
131 // Call implementation
132 //
133 IndexType info = sv_generic(A, piv, B);
134
135 # ifdef CHECK_CXXLAPACK
136 //
137 // Compare results
138 //
139 typename GeMatrix<MA>::NoView A_generic = A;
140 typename DenseVector<VP>::NoView piv_generic = piv;
141 typename GeMatrix<MB>::NoView B_generic = B;
142
143 A = A_org;
144 piv = piv_org;
145 B = B_org;
146
147 IndexType _info = sv_native(A, piv, B);
148
149 bool failed = false;
150 if (! isIdentical(A_generic, A, "A_generic", "A")) {
151 std::cerr << "A_org = " << A_org << std::endl;
152 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
153 std::cerr << "F77LAPACK: A = " << A << std::endl;
154 failed = true;
155 }
156
157 if (! isIdentical(piv_generic, piv, "piv_generic", "piv")) {
158 std::cerr << "piv_org = " << piv_org << std::endl;
159 std::cerr << "CXXLAPACK: piv_generic = " << piv_generic << std::endl;
160 std::cerr << "F77LAPACK: piv = " << piv << std::endl;
161 failed = true;
162 }
163
164 if (! isIdentical(B_generic, B, "B_generic", "B")) {
165 std::cerr << "B_org = " << B_org << std::endl;
166 std::cerr << "CXXLAPACK: B_generic = " << B_generic << std::endl;
167 std::cerr << "F77LAPACK: B = " << B << std::endl;
168 failed = true;
169 }
170
171 if (! isIdentical(info, _info, " info", "_info")) {
172 std::cerr << "CXXLAPACK: info = " << info << std::endl;
173 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
174 failed = true;
175 }
176 if (failed) {
177 ASSERT(0);
178 }
179
180 # endif
181
182 return info;
183 }
184
185 template <typename MA, typename VP, typename VB>
186 typename GeMatrix<MA>::IndexType
187 sv(GeMatrix<MA> &A, DenseVector<VP> &piv, DenseVector<VB> &b)
188 {
189 typedef typename DenseVector<VB>::ElementType ElementType;
190 typedef typename DenseVector<VB>::IndexType IndexType;
191
192 const IndexType n = b.length();
193 const StorageOrder order = GeMatrix<MA>::Engine::order;
194
195 GeMatrix<FullStorageView<ElementType, order> > B(n, 1, b, n);
196
197 return sv(A, piv, B);
198 }
199
200 //-- forwarding ----------------------------------------------------------------
201 template <typename MA, typename VP, typename MB>
202 typename MA::IndexType
203 sv(MA &&A, VP &&piv, MB &&B)
204 {
205 typedef typename MA::IndexType IndexType;
206 CHECKPOINT_ENTER;
207 const IndexType info = sv(A, piv, B);
208 CHECKPOINT_LEAVE;
209 return info;
210 }
211
212 } } // namespace lapack, flens
213
214 #endif // FLENS_LAPACK_GESV_SV_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 DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
36 *
37 * -- LAPACK driver 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_SV_TCC
44 #define FLENS_LAPACK_GESV_SV_TCC 1
45
46 #include <flens/blas/blas.h>
47 #include <flens/lapack/lapack.h>
48
49 namespace flens { namespace lapack {
50
51 //== generic lapack implementation =============================================
52 template <typename MA, typename VP, typename MB>
53 typename GeMatrix<MA>::IndexType
54 sv_generic(GeMatrix<MA> &A, DenseVector<VP> &piv, GeMatrix<MB> &B)
55 {
56 typedef typename GeMatrix<MA>::IndexType IndexType;
57
58 IndexType info = 0;
59
60 //
61 // Compute the LU factorization of A.
62 //
63 info = trf(A, piv);
64 if (info==0) {
65 //
66 // Solve the system A*X = B, overwriting B with X.
67 //
68 trs(NoTrans, A, piv, B);
69 }
70 return info;
71 }
72
73 //== interface for native lapack ===============================================
74
75 #ifdef CHECK_CXXLAPACK
76 template <typename MA, typename VP, typename MB>
77 typename GeMatrix<MA>::IndexType
78 sv_native(GeMatrix<MA> &A, DenseVector<VP> &piv, GeMatrix<MB> &B)
79 {
80 typedef typename GeMatrix<MA>::ElementType ElementType;
81 const INTEGER N = A.numRows();
82 const INTEGER NRHS = B.numCols();
83 const INTEGER LDA = A.leadingDimension();
84 const INTEGER LDB = B.leadingDimension();
85 INTEGER INFO;
86
87 if (IsSame<ElementType, DOUBLE>::value) {
88 LAPACK_IMPL(dgesv)(&N,
89 &NRHS,
90 A.data(),
91 &LDA,
92 piv.data(),
93 B.data(),
94 &LDB,
95 &INFO);
96 } else {
97 ASSERT(0);
98 }
99 return INFO;
100 }
101
102 #endif // CHECK_CXXLAPACK
103
104 //== public interface ==========================================================
105
106 template <typename MA, typename VP, typename MB>
107 typename GeMatrix<MA>::IndexType
108 sv(GeMatrix<MA> &A, DenseVector<VP> &piv, GeMatrix<MB> &B)
109 {
110 typedef typename GeMatrix<MA>::IndexType IndexType;
111 //
112 // Test the input parameters
113 //
114 # ifndef NDEBUG
115 ASSERT(A.firstRow()==1);
116 ASSERT(A.firstCol()==1);
117 ASSERT(A.numRows()==A.numCols());
118 ASSERT((piv.inc()>0 && piv.firstIndex()==1)
119 || (piv.inc()<0 && piv.firstIndex()==A.numRows()));
120 ASSERT(B.firstRow()==1);
121 ASSERT(B.firstCol()==1);
122 ASSERT(B.numRows()==A.numRows());
123 # endif
124 //
125 // Make copies of output arguments
126 //
127 typename GeMatrix<MA>::NoView A_org = A;
128 typename DenseVector<VP>::NoView piv_org = piv;
129 typename GeMatrix<MB>::NoView B_org = B;
130 //
131 // Call implementation
132 //
133 IndexType info = sv_generic(A, piv, B);
134
135 # ifdef CHECK_CXXLAPACK
136 //
137 // Compare results
138 //
139 typename GeMatrix<MA>::NoView A_generic = A;
140 typename DenseVector<VP>::NoView piv_generic = piv;
141 typename GeMatrix<MB>::NoView B_generic = B;
142
143 A = A_org;
144 piv = piv_org;
145 B = B_org;
146
147 IndexType _info = sv_native(A, piv, B);
148
149 bool failed = false;
150 if (! isIdentical(A_generic, A, "A_generic", "A")) {
151 std::cerr << "A_org = " << A_org << std::endl;
152 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
153 std::cerr << "F77LAPACK: A = " << A << std::endl;
154 failed = true;
155 }
156
157 if (! isIdentical(piv_generic, piv, "piv_generic", "piv")) {
158 std::cerr << "piv_org = " << piv_org << std::endl;
159 std::cerr << "CXXLAPACK: piv_generic = " << piv_generic << std::endl;
160 std::cerr << "F77LAPACK: piv = " << piv << std::endl;
161 failed = true;
162 }
163
164 if (! isIdentical(B_generic, B, "B_generic", "B")) {
165 std::cerr << "B_org = " << B_org << std::endl;
166 std::cerr << "CXXLAPACK: B_generic = " << B_generic << std::endl;
167 std::cerr << "F77LAPACK: B = " << B << std::endl;
168 failed = true;
169 }
170
171 if (! isIdentical(info, _info, " info", "_info")) {
172 std::cerr << "CXXLAPACK: info = " << info << std::endl;
173 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
174 failed = true;
175 }
176 if (failed) {
177 ASSERT(0);
178 }
179
180 # endif
181
182 return info;
183 }
184
185 template <typename MA, typename VP, typename VB>
186 typename GeMatrix<MA>::IndexType
187 sv(GeMatrix<MA> &A, DenseVector<VP> &piv, DenseVector<VB> &b)
188 {
189 typedef typename DenseVector<VB>::ElementType ElementType;
190 typedef typename DenseVector<VB>::IndexType IndexType;
191
192 const IndexType n = b.length();
193 const StorageOrder order = GeMatrix<MA>::Engine::order;
194
195 GeMatrix<FullStorageView<ElementType, order> > B(n, 1, b, n);
196
197 return sv(A, piv, B);
198 }
199
200 //-- forwarding ----------------------------------------------------------------
201 template <typename MA, typename VP, typename MB>
202 typename MA::IndexType
203 sv(MA &&A, VP &&piv, MB &&B)
204 {
205 typedef typename MA::IndexType IndexType;
206 CHECKPOINT_ENTER;
207 const IndexType info = sv(A, piv, B);
208 CHECKPOINT_LEAVE;
209 return info;
210 }
211
212 } } // namespace lapack, flens
213
214 #endif // FLENS_LAPACK_GESV_SV_TCC