30char matrice_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Matrice/matrice.C,v 1.19 2014/10/13 08:53:07 j_novak Exp $" ;
151#include "proto_f77.h"
157 if (
std != 0x0)
delete std ;
166 if (
lu != 0x0)
delete lu ;
221 if (source.
lu != 0x0)
lu =
new Tbl(*source.
lu) ;
235 if (source.
get_etat() == ETATZERO) {
239 assert( source.
get_etat() == ETATQCQ ) ;
241 for (
int i=0; i<n; i++) {
242 std->
t[i] = source.
t[i] ;
278 switch (source.
etat) {
289 if (source.
std != 0x0)
292 if (source.
band != 0x0) {
298 if (source.
lu != 0x0) {
311 switch (source.
etat) {
322 assert (source.
t != 0x0) ;
330ostream& operator<< (ostream& flux,
const Matrice & source) {
333 flux <<
"Null matrix. " << endl ;
336 flux <<
"Undefined matrix. " << endl ;
341 flux <<
"Matrix " << nbl <<
" * " << nbc << endl ;
342 for (
int i=0 ; i<nbl ; i++) {
343 for (
int j=0 ; j<nbc ; j++)
344 flux << (*source.
std)(i, j) <<
" " ;
352 flux <<
"Matrix : " << source.
ku <<
" upper diags. and "
353 << source.
kl <<
" lower diags." << endl ;
357 if ((source.
lu != 0x0) && (source.
lu->
get_etat() != ETATNONDEF))
358 flux <<
"LU factorization done." << endl ;
365 if (
band != 0x0) return ;
376 for (
int i=0 ; i<u ; i++)
377 for (
int j=u-i ; j<n ; j++)
378 band->
set(j*ldab+i+l) = (*this)(j-u+i, j) ;
380 for (
int j=0 ; j<n ; j++)
381 band->
set(j*ldab+u+l) = (*this)(j, j) ;
383 for (
int i=u+1 ; i<u+l+1 ; i++)
384 for (
int j=0 ; j<n-i+u ; j++)
385 band->
set(j*ldab+i+l) = (*this) (i+j-u, j) ;
417 F77_dgetrf(&n, &n,
lu->
t, &ldab,
permute->
t, &info) ;
443 F77_dgbtrs(trans, &n, &
kl, &
ku, &nrhs,
lu->
t,
449 F77_dgetrs(trans, &n, &nrhs,
lu->
t, &ldab,
permute->
t,
450 res.
t, &ldb, &info) ;
459 assert (
etat != ETATNONDEF) ;
460 assert (
std != 0x0) ;
462 const char* jobvl =
"N" ;
463 const char* jobvr =
"N" ;
468 double* a =
new double [n*n] ;
469 for (
int i=0 ; i<n*n ; i++)
473 double* wr =
new double[n] ;
474 double* wi =
new double[n] ;
482 double* work =
new double[ldwork] ;
486 F77_dgeev(jobvl, jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
487 work, &ldwork, &info) ;
492 for (
int i=0 ; i<n ; i++) {
493 result.
set(0, i) = wr[n-i-1] ;
494 result.
set(1, i) = wi[n-i-1] ;
509 assert (
etat != ETATNONDEF) ;
510 assert (
std != 0x0) ;
512 const char* jobvl =
"V" ;
513 const char* jobvr =
"N" ;
518 double* a =
new double [n*n] ;
519 for (
int i=0 ; i<n*n ; i++)
523 double* wr =
new double[n] ;
524 double* wi =
new double[n] ;
527 double* vl =
new double[ldvl*ldvl] ;
532 double* work =
new double[ldwork] ;
536 F77_dgeev(jobvl, jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
537 work, &ldwork, &info) ;
544 for (
int i=0 ; i<n ; i++)
545 for (
int j=0 ; j<n ; j++) {
546 res.
set(j,n-i-1) = vl[conte] ;
567 for (
int i = 0 ; i<n ; i++)
569 result *= valp(0, i) ;
571 result*= valp(0, i)*valp(0, i)+valp(1, i)*valp(1, i) ;
585 if (
etat == ETATZERO) {
589 assert(
etat == ETATQCQ) ;
591 for (
int i=0; i<nbc; i++) {
592 for (
int j=0; j<nbl; j++) {
593 resu.
set(i,j) = (*std)(j,i) ;
603 assert((
std != 0x0)&&(a.
std != 0x0)) ;
608 assert((
std != 0x0)&&(a.
std != 0x0)) ;
635 assert((a.
std != 0x0) && (b.
std != 0x0)) ;
641 assert((a.
std != 0x0) && (b.
std != 0x0)) ;
647 assert(a.
std != 0x0) ;
653 assert(a.
std != 0x0) ;
667 assert( nbca == nblb ) ;
675 assert( aa.
get_etat() == ETATQCQ ) ;
676 assert( bb.
get_etat() == ETATQCQ ) ;
678 for (
int i=0; i<nbla; i++) {
679 for (
int j=0; j<nbcb; j++) {
681 for (
int k=0; k<nbca; k++) {
682 sum += aa(i,k) * bb(k, j) ;
684 resu.
set(i,j) = sum ;
695 assert(a.
std != 0x0) ;
Basic integer array class.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
int get_etat() const
Gives the logical state.
int * t
The array of int 's.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void operator=(double x)
Sets all the element of *std to x.
void operator+=(const Matrice &)
Addition of a Matrice to this.
Tbl val_propre() const
Returns the eigenvalues of the matrix, calculated using LAPACK.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
int ku
Number of upper-diagonals in the band representation.
double & set(int j, int i)
Read/write of a particuliar element.
int get_etat() const
Returns the logical state.
int kl
Number of lower-diagonals in the band representation.
void operator/=(double)
Division of this by a double.
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined state).
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
Itbl * permute
Pointer on the second array of the LU-representation.
Tbl * band
Pointer on the array of the band representation of a square matrix.
Tbl * lu
Pointer on the first array of the LU-representation.
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
void del_t()
Logical destructor : dellocates the memory of the various used representations.
double determinant() const
Computes the determinant of the matrix, using LAPACK and the standard decomposition.
int get_dim(int i) const
Returns the dimension of the matrix.
Tbl * std
Pointer on the array of the standard representation.
Matrice transpose() const
Computes the transpose matrix.
void del_deriv()
Deletes the (mutable) derived members: band, lu, permute.
void set_band(int up, int low) const
Calculate the band storage of *std.
void operator*=(double)
Multiplication of this by a double.
int etat
logical state (ETATZERO, ETATQCQ or ETATNONDEF)
void operator-=(const Matrice &)
Subtraction of a Matrice to this.
Matrice(int size1, int size2)
Standard constructor.
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Matrice vect_propre() const
Returns the eigenvectors of the matrix, calculated using LAPACK.
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
int get_etat() const
Gives the logical state.
int etat
logical state (ETATNONDEF, ETATQCQ or ETATZERO).
void annule_hard()
Sets the Tbl to zero in a hard way.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int i)
Read/write of a particular element (index i) (1D case)
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined).
int get_taille() const
Gives the total size (ie dim.taille)
double * t
The array of double.
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
Cmp operator-(const Cmp &)
- Cmp
Cmp operator/(const Cmp &, const Cmp &)
Cmp / Cmp.
Cmp operator+(const Cmp &)