154 #include "proto_f77.h" 160 if (
std != 0x0)
delete std ;
169 if (
lu != 0x0)
delete lu ;
179 std->set_etat_qcq() ;
185 std->set_etat_zero() ;
191 if (
std != 0x0)
std->set_etat_nondef() ;
197 std->set_etat_qcq() ;
201 for (
int i=0 ; i<
std->get_taille() ; i++)
224 if (source.
lu != 0x0)
lu =
new Tbl(*source.
lu) ;
238 if (source.
get_etat() == ETATZERO) {
239 std->set_etat_zero() ;
242 assert( source.
get_etat() == ETATQCQ ) ;
243 std->set_etat_qcq() ;
244 for (
int i=0; i<n; i++) {
245 std->t[i] = source.
t[i] ;
264 return std->get_dim(i) ;
281 switch (source.
etat) {
292 if (source.
std != 0x0)
295 if (source.
band != 0x0) {
301 if (source.
lu != 0x0) {
311 assert (
std->get_dim(0) == source.
get_dim(0)) ;
312 assert (
std->get_dim(1) == source.
get_dim(1)) ;
314 switch (source.
etat) {
325 assert (source.
t != 0x0) ;
333 ostream& operator<< (ostream& flux,
const Matrice & source) {
336 flux <<
"Null matrix. " << endl ;
339 flux <<
"Undefined matrix. " << endl ;
344 flux <<
"Matrix " << nbl <<
" * " << nbc << endl ;
345 for (
int i=0 ; i<nbl ; i++) {
346 for (
int j=0 ; j<nbc ; j++)
347 flux << (*source.
std)(i, j) <<
" " ;
355 flux <<
"Matrix : " << source.
ku <<
" upper diags. and " 356 << source.
kl <<
" lower diags." << endl ;
360 if ((source.
lu != 0x0) && (source.
lu->
get_etat() != ETATNONDEF))
361 flux <<
"LU factorization done." << endl ;
368 if (
band != 0x0) return ;
370 int n =
std->get_dim(0) ;
371 assert (n ==
std->get_dim(1)) ;
379 for (
int i=0 ; i<u ; i++)
380 for (
int j=u-i ; j<n ; j++)
381 band->
set(j*ldab+i+l) = (*this)(j-u+i, j) ;
383 for (
int j=0 ; j<n ; j++)
384 band->
set(j*ldab+u+l) = (*this)(j, j) ;
386 for (
int i=u+1 ; i<u+l+1 ; i++)
387 for (
int j=0 ; j<n-i+u ; j++)
388 band->
set(j*ldab+i+l) = (*this) (i+j-u, j) ;
402 int n =
std->get_dim(0) ;
416 assert (
std->get_etat() == ETATQCQ) ;
420 F77_dgetrf(&n, &n,
lu->
t, &ldab,
permute->
t, &info) ;
446 F77_dgbtrs(trans, &n, &
kl, &
ku, &nrhs,
lu->
t,
452 F77_dgetrs(trans, &n, &nrhs,
lu->
t, &ldab,
permute->
t,
453 res.
t, &ldb, &info) ;
462 assert (
etat != ETATNONDEF) ;
463 assert (
std != 0x0) ;
465 const char* jobvl =
"N" ;
466 const char* jobvr =
"N" ;
471 double* a =
new double [n*n] ;
472 for (
int i=0 ; i<n*n ; i++)
476 double* wr =
new double[n] ;
477 double* wi =
new double[n] ;
485 double* work =
new double[ldwork] ;
489 F77_dgeev(jobvl, jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
490 work, &ldwork, &info) ;
495 for (
int i=0 ; i<n ; i++) {
496 result.
set(0, i) = wr[n-i-1] ;
497 result.
set(1, i) = wi[n-i-1] ;
512 assert (
etat != ETATNONDEF) ;
513 assert (
std != 0x0) ;
515 const char* jobvl =
"V" ;
516 const char* jobvr =
"N" ;
521 double* a =
new double [n*n] ;
522 for (
int i=0 ; i<n*n ; i++)
526 double* wr =
new double[n] ;
527 double* wi =
new double[n] ;
530 double* vl =
new double[ldvl*ldvl] ;
535 double* work =
new double[ldwork] ;
539 F77_dgeev(jobvl, jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr,
540 work, &ldwork, &info) ;
547 for (
int i=0 ; i<n ; i++)
548 for (
int j=0 ; j<n ; j++) {
549 res.
set(j,n-i-1) = vl[conte] ;
570 for (
int i = 0 ; i<n ; i++)
572 result *= valp(0, i) ;
574 result*= valp(0, i)*valp(0, i)+valp(1, i)*valp(1, i) ;
583 int nbl =
std->get_dim(1) ;
584 int nbc =
std->get_dim(0) ;
588 if (
etat == ETATZERO) {
592 assert(
etat == ETATQCQ) ;
594 for (
int i=0; i<nbc; i++) {
595 for (
int j=0; j<nbl; j++) {
596 resu.
set(i,j) = (*std)(j,i) ;
606 assert((
std != 0x0)&&(a.
std != 0x0)) ;
611 assert((
std != 0x0)&&(a.
std != 0x0)) ;
638 assert((a.
std != 0x0) && (b.
std != 0x0)) ;
644 assert((a.
std != 0x0) && (b.
std != 0x0)) ;
650 assert(a.
std != 0x0) ;
656 assert(a.
std != 0x0) ;
670 assert( nbca == nblb ) ;
678 assert( aa.
get_etat() == ETATQCQ ) ;
679 assert( bb.
get_etat() == ETATQCQ ) ;
681 for (
int i=0; i<nbla; i++) {
682 for (
int j=0; j<nbcb; j++) {
684 for (
int k=0; k<nbca; k++) {
685 sum += aa(i,k) * bb(k, j) ;
687 resu.
set(i,j) = sum ;
698 assert(a.
std != 0x0) ;
Tbl * lu
Pointer on the first array of the LU-representation.
void operator/=(double)
Division of this by a double.
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Base_val operator*(const Base_val &, const Base_val &)
This operator is used when calling multiplication or division of Valeur .
int get_etat() const
Returns the logical state.
Basic integer array class.
int get_etat() const
Gives the logical state.
Cmp operator/(const Cmp &, const Cmp &)
Cmp / Cmp.
Tbl * band
Pointer on the array of the band representation of a square matrix.
void operator-=(const Matrice &)
Subtraction of a Matrice to this.
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined state).
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
int ku
Number of upper-diagonals in the band representation.
int get_ndim() const
Gives the number of dimensions (ie dim.ndim)
void operator+=(const Matrice &)
Addition of a Matrice to this.
Tbl val_propre() const
Returns the eigenvalues of the matrix, calculated using LAPACK.
Cmp operator+(const Cmp &)
double determinant() const
Computes the determinant of the matrix, using LAPACK and the standard decomposition.
int kl
Number of lower-diagonals in the band representation.
double * t
The array of double.
void operator*=(double)
Multiplication of this by a double.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Itbl * permute
Pointer on the second array of the LU-representation.
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
int get_dim(int i) const
Gives the i-th dimension (ie dim.dim[i])
int get_dim(int i) const
Returns the dimension of the matrix.
Matrice transpose() const
Computes the transpose matrix.
int etat
logical state (ETATNONDEF, ETATQCQ or ETATZERO).
double & set(int j, int i)
Read/write of a particuliar element.
Tbl * std
Pointer on the array of the standard representation.
int get_etat() const
Gives the logical state.
void set_band(int up, int low) const
Calculate the band storage of *std.
void del_deriv()
Deletes the (mutable) derived members: band, lu, permute.
int get_taille() const
Gives the total size (ie dim.taille)
void operator=(double x)
Sets all the element of *std to x.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
int etat
logical state (ETATZERO, ETATQCQ or ETATNONDEF)
Cmp operator-(const Cmp &)
- Cmp
int * t
The array of int 's.
void annule_hard()
Sets the Tbl to zero in a hard way.
Matrice(int size1, int size2)
Standard constructor.
Matrice vect_propre() const
Returns the eigenvectors of the matrix, calculated using LAPACK.
void del_t()
Logical destructor : dellocates the memory of the various used representations.