BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
fn_geev.hpp
Go to the documentation of this file.
1/*
2 * BELFEM -- The Berkeley Lab Finite Element Framework
3 * Copyright (c) 2026, The Regents of the University of California, through
4 * Lawrence Berkeley National Laboratory (subject to receipt of any required
5 * approvals from the U.S. Dept. of Energy). All rights reserved.
6 *
7 * Developers: Christian Messe, Gregory Giard
8 *
9 * See the top-level LICENSE file for the complete license and disclaimer.
10 */
11
21
22#ifndef BELFEM_FN_GEEV_HPP
23#define BELFEM_FN_GEEV_HPP
24
25#include "assert.hpp"
26#include "lapacktools.hpp"
27#include "cl_Vector.hpp"
28
29//------------------------------------------------------------------------------
30namespace belfem
31{
32 namespace lapack
33 {
34//------------------------------------------------------------------------------
35#ifdef __cplusplus
36 extern "C"
37 {
38#endif
39//------------------------------------------------------------------------------
40
41 // geev particularities: the real flavors return the eigenvalues as
42 // separate real and imaginary arrays ( wr, wi ) and have no rwork;
43 // the complex flavors take one complex w plus a REAL rwork of size
44 // 2*n. The two job chars carry two hidden Fortran length arguments.
45
46 void
47 sgeev_( char * jobvl,
48 char * jobvr,
49 int_t * n,
50 float * a,
51 int_t * lda,
52 float * wr,
53 float * wi,
54 float * vl,
55 int_t * ldvl,
56 float * vr,
57 int_t * ldvr,
58 float * work,
59 int_t * lwork,
60 int_t * info,
61 fortran_charlen_t lvl,
62 fortran_charlen_t lvr );
63
64//------------------------------------------------------------------------------
65
66 void
67 dgeev_( char * jobvl,
68 char * jobvr,
69 int_t * n,
70 double * a,
71 int_t * lda,
72 double * wr,
73 double * wi,
74 double * vl,
75 int_t * ldvl,
76 double * vr,
77 int_t * ldvr,
78 double * work,
79 int_t * lwork,
80 int_t * info,
81 fortran_charlen_t lvl,
82 fortran_charlen_t lvr );
83
84//------------------------------------------------------------------------------
85
86 void
87 cgeev_( char * jobvl,
88 char * jobvr,
89 int_t * n,
90 cplx_float_t * a,
91 int_t * lda,
92 cplx_float_t * w,
93 cplx_float_t * vl,
94 int_t * ldvl,
95 cplx_float_t * vr,
96 int_t * ldvr,
97 cplx_float_t * work,
98 int_t * lwork,
99 float * rwork,
100 int_t * info,
101 fortran_charlen_t lvl,
102 fortran_charlen_t lvr );
103
104//------------------------------------------------------------------------------
105
106 void
107 zgeev_( char * jobvl,
108 char * jobvr,
109 int_t * n,
110 cplx_double_t * a,
111 int_t * lda,
112 cplx_double_t * w,
113 cplx_double_t * vl,
114 int_t * ldvl,
115 cplx_double_t * vr,
116 int_t * ldvr,
117 cplx_double_t * work,
118 int_t * lwork,
119 double * rwork,
120 int_t * info,
121 fortran_charlen_t lvl,
122 fortran_charlen_t lvr );
123
124//------------------------------------------------------------------------------
125#ifdef __cplusplus
126 }
127#endif
128//------------------------------------------------------------------------------
129
130 // unified dispatch. The real scratch rwork must hold 2*n entries in
131 // ALL flavors: the real specializations use it as wr/wi and pack the
132 // eigenvalues into the complex w afterwards, the complex ones pass
133 // it through as LAPACK's actual rwork.
134
135 template< typename T >
136 void geev(
137 const char * jobvl,
138 const char * jobvr,
139 const int_t * n,
140 T * a,
141 const int_t * lda,
142 cplx_t< T > * w,
143 T * vl,
144 const int_t * ldvl,
145 T * vr,
146 const int_t * ldvr,
147 T * work,
148 const int_t * lwork,
149 real_t< T > * rwork,
150 int_t * info )
151 {
152 static_assert( dependent_false< T >,
153 "geev not implemented for selected data type" );
154 }
155
156// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157
158 template<>
159 inline void geev(
160 const char * jobvl,
161 const char * jobvr,
162 const int_t * n,
163 float * a,
164 const int_t * lda,
165 std::complex< float > * w,
166 float * vl,
167 const int_t * ldvl,
168 float * vr,
169 const int_t * ldvr,
170 float * work,
171 const int_t * lwork,
172 float * rwork,
173 int_t * info )
174 {
175 sgeev_(
176 const_cast< char * >( jobvl ),
177 const_cast< char * >( jobvr ),
178 const_cast< int_t * >( n ),
179 a,
180 const_cast< int_t * >( lda ),
181 rwork,
182 rwork + *n,
183 vl,
184 const_cast< int_t * >( ldvl ),
185 vr,
186 const_cast< int_t * >( ldvr ),
187 work,
188 const_cast< int_t * >( lwork ),
189 info,
190 1, 1 );
191
192 // pack wr/wi into the complex eigenvalue vector; a workspace
193 // query does not touch wr/wi, so there is nothing to pack
194 if ( *lwork != -1 )
195 {
196 for ( int_t k = 0; k < *n; ++k )
197 {
198 w[ k ] = std::complex< float >( rwork[ k ], rwork[ k + *n ] );
199 }
200 }
201 }
202
203// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204
205 template<>
206 inline void geev(
207 const char * jobvl,
208 const char * jobvr,
209 const int_t * n,
210 double * a,
211 const int_t * lda,
212 std::complex< double > * w,
213 double * vl,
214 const int_t * ldvl,
215 double * vr,
216 const int_t * ldvr,
217 double * work,
218 const int_t * lwork,
219 double * rwork,
220 int_t * info )
221 {
222 dgeev_(
223 const_cast< char * >( jobvl ),
224 const_cast< char * >( jobvr ),
225 const_cast< int_t * >( n ),
226 a,
227 const_cast< int_t * >( lda ),
228 rwork,
229 rwork + *n,
230 vl,
231 const_cast< int_t * >( ldvl ),
232 vr,
233 const_cast< int_t * >( ldvr ),
234 work,
235 const_cast< int_t * >( lwork ),
236 info,
237 1, 1 );
238
239 if ( *lwork != -1 )
240 {
241 for ( int_t k = 0; k < *n; ++k )
242 {
243 w[ k ] = std::complex< double >( rwork[ k ], rwork[ k + *n ] );
244 }
245 }
246 }
247
248// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249
250 template<>
251 inline void geev(
252 const char * jobvl,
253 const char * jobvr,
254 const int_t * n,
255 std::complex< float > * a,
256 const int_t * lda,
257 std::complex< float > * w,
258 std::complex< float > * vl,
259 const int_t * ldvl,
260 std::complex< float > * vr,
261 const int_t * ldvr,
262 std::complex< float > * work,
263 const int_t * lwork,
264 float * rwork,
265 int_t * info )
266 {
267 cgeev_(
268 const_cast< char * >( jobvl ),
269 const_cast< char * >( jobvr ),
270 const_cast< int_t * >( n ),
271 reinterpret_cast< cplx_float_t * >( a ),
272 const_cast< int_t * >( lda ),
273 reinterpret_cast< cplx_float_t * >( w ),
274 reinterpret_cast< cplx_float_t * >( vl ),
275 const_cast< int_t * >( ldvl ),
276 reinterpret_cast< cplx_float_t * >( vr ),
277 const_cast< int_t * >( ldvr ),
278 reinterpret_cast< cplx_float_t * >( work ),
279 const_cast< int_t * >( lwork ),
280 rwork,
281 info,
282 1, 1 );
283 }
284
285// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286
287 template<>
288 inline void geev(
289 const char * jobvl,
290 const char * jobvr,
291 const int_t * n,
292 std::complex< double > * a,
293 const int_t * lda,
294 std::complex< double > * w,
295 std::complex< double > * vl,
296 const int_t * ldvl,
297 std::complex< double > * vr,
298 const int_t * ldvr,
299 std::complex< double > * work,
300 const int_t * lwork,
301 double * rwork,
302 int_t * info )
303 {
304 zgeev_(
305 const_cast< char * >( jobvl ),
306 const_cast< char * >( jobvr ),
307 const_cast< int_t * >( n ),
308 reinterpret_cast< cplx_double_t * >( a ),
309 const_cast< int_t * >( lda ),
310 reinterpret_cast< cplx_double_t * >( w ),
311 reinterpret_cast< cplx_double_t * >( vl ),
312 const_cast< int_t * >( ldvl ),
313 reinterpret_cast< cplx_double_t * >( vr ),
314 const_cast< int_t * >( ldvr ),
315 reinterpret_cast< cplx_double_t * >( work ),
316 const_cast< int_t * >( lwork ),
317 rwork,
318 info,
319 1, 1 );
320 }
321
322//------------------------------------------------------------------------------
323 } /* end namespace lapack */
324
352 template< typename T >
353 int_t
354 geev(
355 Matrix< T > & A,
356 Vector< lapack::cplx_t< T > > & W,
357 Matrix< T > & VL,
358 Matrix< T > & VR,
359 Vector< lapack::real_t< T > > & Work,
360 const char jobvl = 'V',
361 const char jobvr = 'V',
362 const bool AbortOnError = true )
363 {
364 BELFEM_ASSERT( A.n_rows() == A.n_cols(),
365 "Matrix A must be square ( is %lu x %lu )",
366 ( long unsigned int ) A.n_rows(),
367 ( long unsigned int ) A.n_cols() );
368 BELFEM_ASSERT( jobvl == 'V' || jobvl == 'N',
369 "unsupported jobvl flag '%c'", jobvl );
370 BELFEM_ASSERT( jobvr == 'V' || jobvr == 'N',
371 "unsupported jobvr flag '%c'", jobvr );
372
373 // reals per LAPACK work entry
374 constexpr int_t tRealsPerT =
375 std::is_same< T, lapack::real_t< T > >::value ? 1 : 2 ;
376
377 int_t n = ( int_t ) A.n_rows();
378 int_t lda = lapack::leading_dimension( A );
379
380 W.set_size( n );
381
382 // lapack requires ldvl/ldvr >= 1 even when not referenced
383 int_t ldvl = 1 ;
384 if ( jobvl == 'V' )
385 {
386 VL.set_size( n, n );
387 ldvl = lapack::leading_dimension( VL );
388 }
389
390 int_t ldvr = 1 ;
391 if ( jobvr == 'V' )
392 {
393 VR.set_size( n, n );
394 ldvr = lapack::leading_dimension( VR );
395 }
396
397 // minimum work sizes: real flavors need 4n with eigenvectors and
398 // 3n without, the complex ones need 2n
399 int_t lwork = std::max< int_t >( 1, tRealsPerT == 1 ?
400 ( ( jobvl == 'V' || jobvr == 'V' ) ? 4 * n : 3 * n ) : 2 * n );
401
402 int_t info = 0 ;
403
404 // required buffer: lwork entries of T in the head plus 2*n reals
405 if ( static_cast< int_t >( Work.length() ) < tRealsPerT * lwork + 2 * n )
406 {
407 // ask lapack for the optimal size; one T entry plus the tail
408 Work.set_size( tRealsPerT + 2 * n );
409
410 int_t query = -1 ;
411
412 lapack::geev( &jobvl, &jobvr, &n, A.data(), &lda, W.data(),
413 VL.data(), &ldvl, VR.data(), &ldvr,
414 reinterpret_cast< T * >( Work.data() ), &query,
415 Work.data() + tRealsPerT, &info );
416
417 BELFEM_ERROR( info == 0 || ! AbortOnError,
418 "LAPACK geev workspace query has thrown an error: %i", ( int ) info );
419
420 if ( info != 0 ) return info ;
421
422 // the optimal size sits in the real part of the first entry
423 lwork = lapack::work_size( Work( 0 ) );
424
425 Work.set_size( tRealsPerT * lwork + 2 * n );
426 }
427 else
428 {
429 // use the full buffer the caller has provided
430 lwork = ( ( int_t ) Work.length() - 2 * n ) / tRealsPerT ;
431 }
432
433 // work segment in the head, 2*n real scratch in the tail
434 lapack::geev( &jobvl, &jobvr, &n, A.data(), &lda, W.data(),
435 VL.data(), &ldvl, VR.data(), &ldvr,
436 reinterpret_cast< T * >( Work.data() ), &lwork,
437 Work.data() + tRealsPerT * lwork, &info );
438
439 // info > 0 : the QR iteration failed; eigenvalues info+1 .. n
440 // have converged, no eigenvectors were computed
441 BELFEM_ERROR( info == 0 || ! AbortOnError,
442 "LAPACK geev has thrown an error: %i", ( int ) info );
443
444 return info ;
445 }
446
447//------------------------------------------------------------------------------
448
460 template< typename T >
461 int_t
462 geev(
463 Matrix< T > & A,
464 Vector< lapack::cplx_t< T > > & W,
465 Vector< lapack::real_t< T > > & Work,
466 const bool AbortOnError = true )
467 {
468 BELFEM_ASSERT( A.n_rows() == A.n_cols(),
469 "Matrix A must be square ( is %lu x %lu )",
470 ( long unsigned int ) A.n_rows(),
471 ( long unsigned int ) A.n_cols() );
472
473 constexpr int_t tRealsPerT =
474 std::is_same< T, lapack::real_t< T > >::value ? 1 : 2 ;
475
476 char jobv = 'N' ;
477
478 int_t n = ( int_t ) A.n_rows();
479 int_t lda = lapack::leading_dimension( A );
480
481 W.set_size( n );
482
483 // not referenced under jobvl = jobvr = 'N', but ld >= 1 is required
484 int_t ldv = 1 ;
485
486 int_t lwork = std::max< int_t >( 1, tRealsPerT == 1 ? 3 * n : 2 * n );
487
488 int_t info = 0 ;
489
490 if ( static_cast< int_t >( Work.length() ) < tRealsPerT * lwork + 2 * n )
491 {
492 Work.set_size( tRealsPerT + 2 * n );
493
494 int_t query = -1 ;
495
496 lapack::geev( &jobv, &jobv, &n, A.data(), &lda, W.data(),
497 static_cast< T * >( nullptr ), &ldv,
498 static_cast< T * >( nullptr ), &ldv,
499 reinterpret_cast< T * >( Work.data() ), &query,
500 Work.data() + tRealsPerT, &info );
501
502 BELFEM_ERROR( info == 0 || ! AbortOnError,
503 "LAPACK geev workspace query has thrown an error: %i", ( int ) info );
504
505 if ( info != 0 ) return info ;
506
507 lwork = lapack::work_size( Work( 0 ) );
508
509 Work.set_size( tRealsPerT * lwork + 2 * n );
510 }
511 else
512 {
513 lwork = ( ( int_t ) Work.length() - 2 * n ) / tRealsPerT ;
514 }
515
516 lapack::geev( &jobv, &jobv, &n, A.data(), &lda, W.data(),
517 static_cast< T * >( nullptr ), &ldv,
518 static_cast< T * >( nullptr ), &ldv,
519 reinterpret_cast< T * >( Work.data() ), &lwork,
520 Work.data() + tRealsPerT * lwork, &info );
521
522 BELFEM_ERROR( info == 0 || ! AbortOnError,
523 "LAPACK geev has thrown an error: %i", ( int ) info );
524
525 return info ;
526 }
527
528//------------------------------------------------------------------------------
529}
530#endif //BELFEM_FN_GEEV_HPP
BELFEM_ERROR
#define BELFEM_ERROR(aCheck,...)
Definition assert.hpp:264
BELFEM_ASSERT
#define BELFEM_ASSERT(aCheck,...)
Definition assert.hpp:244
cl_Vector.hpp
belfem::Matrix
Dense column-major matrix.
Definition cl_BZ_Matrix.hpp:28
belfem::Matrix::set_size
void set_size(const size_t aNumRows, const size_t aNumCols)
Definition cl_AR_Matrix.hpp:186
belfem::Matrix::n_rows
size_t n_rows() const
Definition cl_AR_Matrix.hpp:205
belfem::Matrix::n_cols
size_t n_cols() const
Definition cl_AR_Matrix.hpp:213
belfem::Matrix::data
T * data()
Definition cl_AR_Matrix.hpp:135
belfem::Vector
Column vector.
Definition cl_BZ_Vector.hpp:41
lapacktools.hpp
Shared helpers for the LAPACK wrappers, such as leading_dimension().
belfem::lapack
Definition fn_gees.hpp:33
belfem::lapack::geev
void geev(const char *jobvl, const char *jobvr, const int_t *n, T *a, const int_t *lda, cplx_t< T > *w, T *vl, const int_t *ldvl, T *vr, const int_t *ldvr, T *work, const int_t *lwork, real_t< T > *rwork, int_t *info)
Definition fn_geev.hpp:136
belfem::lapack::cgeev_
void cgeev_(char *jobvl, char *jobvr, int_t *n, cplx_float_t *a, int_t *lda, cplx_float_t *w, cplx_float_t *vl, int_t *ldvl, cplx_float_t *vr, int_t *ldvr, cplx_float_t *work, int_t *lwork, float *rwork, int_t *info, fortran_charlen_t lvl, fortran_charlen_t lvr)
belfem::lapack::zgeev_
void zgeev_(char *jobvl, char *jobvr, int_t *n, cplx_double_t *a, int_t *lda, cplx_double_t *w, cplx_double_t *vl, int_t *ldvl, cplx_double_t *vr, int_t *ldvr, cplx_double_t *work, int_t *lwork, double *rwork, int_t *info, fortran_charlen_t lvl, fortran_charlen_t lvr)
belfem::lapack::dependent_false
constexpr bool dependent_false
Definition lapacktools.hpp:49
belfem::lapack::cplx_float_t
float cplx_float_t
Definition lapacktools.hpp:90
belfem::lapack::leading_dimension
int_t leading_dimension(const Vector< T > &A)
logical length of a vector operand, as passed to LAPACK as LDB; vector storage is contiguous under bo...
Definition lapacktools.hpp:143
belfem::lapack::work_size
int_t work_size(const float &aValue)
lapack reports the optimal work size in the first entry of the work array, which stays real valued fo...
Definition lapacktools.hpp:191
belfem::lapack::fortran_charlen_t
size_t fortran_charlen_t
Definition lapacktools.hpp:100
belfem::lapack::sgeev_
void sgeev_(char *jobvl, char *jobvr, int_t *n, float *a, int_t *lda, float *wr, float *wi, float *vl, int_t *ldvl, float *vr, int_t *ldvr, float *work, int_t *lwork, int_t *info, fortran_charlen_t lvl, fortran_charlen_t lvr)
belfem::lapack::dgeev_
void dgeev_(char *jobvl, char *jobvr, int_t *n, double *a, int_t *lda, double *wr, double *wi, double *vl, int_t *ldvl, double *vr, int_t *ldvr, double *work, int_t *lwork, int_t *info, fortran_charlen_t lvl, fortran_charlen_t lvr)
belfem::lapack::cplx_double_t
double cplx_double_t
Definition lapacktools.hpp:91
belfem::lapack::real_t
typename real_type< T >::type real_t
Definition lapacktools.hpp:63
belfem::lapack::cplx_t
typename cplx_type< T >::type cplx_t
Definition lapacktools.hpp:75
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::value
std::pair< real, unit > value
Definition typedefs.hpp:74
belfem::geev
int_t geev(Matrix< T > &A, Vector< lapack::cplx_t< T > > &W, Matrix< T > &VL, Matrix< T > &VR, Vector< lapack::real_t< T > > &Work, const char jobvl='V', const char jobvr='V', const bool AbortOnError=true)
eigenvalues and eigenvectors of a general square matrix, A * v = lambda * v, via LAPACK ?...
Definition fn_geev.hpp:354
belfem::int_t
int32_t int_t
Definition typedefs.hpp:51