1 //=================================================================================================
2 /*!
3 // \file blaze/math/lapack/clapack/hesv.h
4 // \brief Header file for the CLAPACK hesv wrapper functions
5 //
6 // Copyright (C) 2012-2020 Klaus Iglberger - All Rights Reserved
7 //
8 // This file is part of the Blaze library. You can redistribute it and/or modify it under
9 // the terms of the New (Revised) BSD License. Redistribution and use in source and binary
10 // forms, with or without modification, are permitted provided that the following conditions
11 // are met:
12 //
13 // 1. Redistributions of source code must retain the above copyright notice, this list of
14 // conditions and the following disclaimer.
15 // 2. Redistributions in binary form must reproduce the above copyright notice, this list
16 // of conditions and the following disclaimer in the documentation and/or other materials
17 // provided with the distribution.
18 // 3. Neither the names of the Blaze development group nor the names of its contributors
19 // may be used to endorse or promote products derived from this software without specific
20 // prior written permission.
21 //
22 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
23 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
24 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
25 // SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
26 // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
27 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
28 // BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30 // ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
31 // DAMAGE.
32 */
33 //=================================================================================================
34
35 #ifndef _BLAZE_MATH_LAPACK_CLAPACK_HESV_H_
36 #define _BLAZE_MATH_LAPACK_CLAPACK_HESV_H_
37
38
39 //*************************************************************************************************
40 // Includes
41 //*************************************************************************************************
42
43 #include <blaze/math/blas/Types.h>
44 #include <blaze/util/Complex.h>
45 #include <blaze/util/StaticAssert.h>
46 #include <blaze/util/Types.h>
47
48
49 //=================================================================================================
50 //
51 // LAPACK FORWARD DECLARATIONS
52 //
53 //=================================================================================================
54
55 //*************************************************************************************************
56 /*! \cond BLAZE_INTERNAL */
57 #if !defined(INTEL_MKL_VERSION)
58 extern "C" {
59
60 void chesv_( char* uplo, blaze::blas_int_t* n, blaze::blas_int_t* nrhs, float* A,
61 blaze::blas_int_t* lda, blaze::blas_int_t* ipiv, float* b, blaze::blas_int_t* ldb,
62 float* work, blaze::blas_int_t* lwork, blaze::blas_int_t* info,
63 blaze::fortran_charlen_t nuplo );
64 void zhesv_( char* uplo, blaze::blas_int_t* n, blaze::blas_int_t* nrhs, double* A,
65 blaze::blas_int_t* lda, blaze::blas_int_t* ipiv, double* b, blaze::blas_int_t* ldb,
66 double* work, blaze::blas_int_t* lwork, blaze::blas_int_t* info,
67 blaze::fortran_charlen_t nuplo );
68
69 }
70 #endif
71 /*! \endcond */
72 //*************************************************************************************************
73
74
75
76
77 namespace blaze {
78
79 //=================================================================================================
80 //
81 // LAPACK HERMITIAN INDEFINITE LINEAR SYSTEM FUNCTIONS (HESV)
82 //
83 //=================================================================================================
84
85 //*************************************************************************************************
86 /*!\name LAPACK Hermitian indefinite linear system functions (hesv) */
87 //@{
88 void hesv( char uplo, blas_int_t n, blas_int_t nrhs, complex<float>* A, blas_int_t lda,
89 blas_int_t* ipiv, complex<float>* B, blas_int_t ldb, complex<float>* work,
90 blas_int_t lwork, blas_int_t* info );
91
92 void hesv( char uplo, blas_int_t n, blas_int_t nrhs, complex<double>* A, blas_int_t lda,
93 blas_int_t* ipiv, complex<double>* B, blas_int_t ldb, complex<double>* work,
94 blas_int_t lwork, blas_int_t* info );
95 //@}
96 //*************************************************************************************************
97
98
99 //*************************************************************************************************
100 /*!\brief LAPACK kernel for solving a Hermitian indefinite single precision complex linear system
101 // of equations (\f$ A*X=B \f$).
102 // \ingroup lapack_solver
103 //
104 // \param uplo \c 'L' to use the lower part of the matrix, \c 'U' to use the upper part.
105 // \param n The number of rows/columns of matrix \a A \f$[0..\infty)\f$.
106 // \param nrhs The number of right-hand side vectors \f$[0..\infty)\f$.
107 // \param A Pointer to the first element of the single precision complex column-major square matrix.
108 // \param lda The total number of elements between two columns of matrix \a A \f$[0..\infty)\f$.
109 // \param ipiv Auxiliary array of size \a n for the pivot indices.
110 // \param B Pointer to the first element of the column-major matrix.
111 // \param ldb The total number of elements between two columns of matrix \a B \f$[0..\infty)\f$.
112 // \param work Auxiliary array; size >= max( 1, \a lwork ).
113 // \param lwork The dimension of the array \a work; size >= max( 1, \a n ).
114 // \param info Return code of the function call.
115 // \return void
116 //
117 // This function uses the LAPACK chesv() function to compute the solution to the Hermitian
118 // indefinite system of linear equations \f$ A*X=B \f$, where \a A is a \a n-by-\a n matrix and
119 // \a X and \a B are \a n-by-\a nrhs matrices.
120 //
121 // The Bunch-Kaufman decomposition is used to factor \a A as
122
123 \f[ A = U D U^{H} \texttt{ (if uplo = 'U'), or }
124 A = L D L^{H} \texttt{ (if uplo = 'L'), } \f]
125
126 // where \c U (or \c L) is a product of permutation and unit upper (lower) triangular matrices,
127 // and \c D is Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The resulting
128 // decomposition is stored within \a A: In case \a uplo is set to \c 'L' the result is stored in
129 // the lower part of the matrix and the upper part remains untouched, in case \a uplo is set to
130 // \c 'U' the result is stored in the upper part and the lower part remains untouched. The factored
131 // form of \a A is then used to solve the system of equations.
132 //
133 // The \a info argument provides feedback on the success of the function call:
134 //
135 // - = 0: The function finished successfully.
136 // - < 0: If info = -i, the i-th argument had an illegal value.
137 // - > 0: If info = i, the decomposition has been completed, but since factor D(i,i) is exactly
138 // zero the solution could not be computed.
139 //
140 // For more information on the chesv() function, see the LAPACK online documentation browser:
141 //
142 // http://www.netlib.org/lapack/explore-html/
143 //
144 // \note This function can only be used if a fitting LAPACK library, which supports this function,
145 // is available and linked to the executable. Otherwise a call to this function will result in a
146 // linker error.
147 */
hesv(char uplo,blas_int_t n,blas_int_t nrhs,complex<float> * A,blas_int_t lda,blas_int_t * ipiv,complex<float> * B,blas_int_t ldb,complex<float> * work,blas_int_t lwork,blas_int_t * info)148 inline void hesv( char uplo, blas_int_t n, blas_int_t nrhs, complex<float>* A, blas_int_t lda,
149 blas_int_t* ipiv, complex<float>* B, blas_int_t ldb, complex<float>* work,
150 blas_int_t lwork, blas_int_t* info )
151 {
152 BLAZE_STATIC_ASSERT( sizeof( complex<float> ) == 2UL*sizeof( float ) );
153
154 #if defined(INTEL_MKL_VERSION)
155 BLAZE_STATIC_ASSERT( sizeof( MKL_INT ) == sizeof( blas_int_t ) );
156 BLAZE_STATIC_ASSERT( sizeof( MKL_Complex8 ) == sizeof( complex<float> ) );
157 using ET = MKL_Complex8;
158 #else
159 using ET = float;
160 #endif
161
162 chesv_( &uplo, &n, &nrhs, reinterpret_cast<ET*>( A ), &lda, ipiv,
163 reinterpret_cast<ET*>( B ), &ldb, reinterpret_cast<ET*>( work ), &lwork, info
164 #if !defined(INTEL_MKL_VERSION)
165 , blaze::fortran_charlen_t(1)
166 #endif
167 );
168 }
169 //*************************************************************************************************
170
171
172 //*************************************************************************************************
173 /*!\brief LAPACK kernel for solving a Hermitian indefinite double precision complex linear system
174 // of equations (\f$ A*X=B \f$).
175 // \ingroup lapack_solver
176 //
177 // \param uplo \c 'L' to use the lower part of the matrix, \c 'U' to use the upper part.
178 // \param n The number of rows/columns of matrix \a A \f$[0..\infty)\f$.
179 // \param nrhs The number of right-hand side vectors \f$[0..\infty)\f$.
180 // \param A Pointer to the first element of the double precision complex column-major square matrix.
181 // \param lda The total number of elements between two columns of matrix \a A \f$[0..\infty)\f$.
182 // \param ipiv Auxiliary array of size \a n for the pivot indices.
183 // \param B Pointer to the first element of the column-major matrix.
184 // \param ldb The total number of elements between two columns of matrix \a B \f$[0..\infty)\f$.
185 // \param work Auxiliary array; size >= max( 1, \a lwork ).
186 // \param lwork The dimension of the array \a work; size >= max( 1, \a n ).
187 // \param info Return code of the function call.
188 // \return void
189 //
190 // This function uses the LAPACK zhesv() function to compute the solution to the Hermitian
191 // indefinite system of linear equations \f$ A*X=B \f$, where \a A is a \a n-by-\a n matrix and
192 // \a X and \a B are \a n-by-\a nrhs matrices.
193 //
194 // The Bunch-Kaufman decomposition is used to factor \a A as
195
196 \f[ A = U D U^{H} \texttt{ (if uplo = 'U'), or }
197 A = L D L^{H} \texttt{ (if uplo = 'L'), } \f]
198
199 // where \c U (or \c L) is a product of permutation and unit upper (lower) triangular matrices,
200 // and \c D is Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The resulting
201 // decomposition is stored within \a A: In case \a uplo is set to \c 'L' the result is stored in
202 // the lower part of the matrix and the upper part remains untouched, in case \a uplo is set to
203 // \c 'U' the result is stored in the upper part and the lower part remains untouched. The factored
204 // form of \a A is then used to solve the system of equations.
205 //
206 // The \a info argument provides feedback on the success of the function call:
207 //
208 // - = 0: The function finished successfully.
209 // - < 0: If info = -i, the i-th argument had an illegal value.
210 // - > 0: If info = i, the decomposition has been completed, but since factor D(i,i) is exactly
211 // zero the solution could not be computed.
212 //
213 // For more information on the zhesv() function, see the LAPACK online documentation browser:
214 //
215 // http://www.netlib.org/lapack/explore-html/
216 //
217 // \note This function can only be used if a fitting LAPACK library, which supports this function,
218 // is available and linked to the executable. Otherwise a call to this function will result in a
219 // linker error.
220 */
hesv(char uplo,blas_int_t n,blas_int_t nrhs,complex<double> * A,blas_int_t lda,blas_int_t * ipiv,complex<double> * B,blas_int_t ldb,complex<double> * work,blas_int_t lwork,blas_int_t * info)221 inline void hesv( char uplo, blas_int_t n, blas_int_t nrhs, complex<double>* A, blas_int_t lda,
222 blas_int_t* ipiv, complex<double>* B, blas_int_t ldb, complex<double>* work,
223 blas_int_t lwork, blas_int_t* info )
224 {
225 BLAZE_STATIC_ASSERT( sizeof( complex<double> ) == 2UL*sizeof( double ) );
226
227 #if defined(INTEL_MKL_VERSION)
228 BLAZE_STATIC_ASSERT( sizeof( MKL_INT ) == sizeof( blas_int_t ) );
229 BLAZE_STATIC_ASSERT( sizeof( MKL_Complex16 ) == sizeof( complex<double> ) );
230 using ET = MKL_Complex16;
231 #else
232 using ET = double;
233 #endif
234
235 zhesv_( &uplo, &n, &nrhs, reinterpret_cast<ET*>( A ), &lda, ipiv,
236 reinterpret_cast<ET*>( B ), &ldb, reinterpret_cast<ET*>( work ), &lwork, info
237 #if !defined(INTEL_MKL_VERSION)
238 , blaze::fortran_charlen_t(1)
239 #endif
240 );
241 }
242 //*************************************************************************************************
243
244 } // namespace blaze
245
246 #endif
247