lanczos.c

cfd求解器使用与gmsh网格的求解

#define RCSID "$Id: Lanczos.c,v 1.34 2006/02/26 00:42:54 geuzaine Exp $"
/*
 * Copyright (C) 1997-2006 P. Dular, C. Geuzaine
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA.
 *
 * Please report all bugs and problems to <[e-mail]>.
 *
 * Contributor(s):
 *   Benoit Meys
 *   Andre Nicolet
 *   [e-mail]
 */

#include "GetDP.h"
#include "DofData.h"
#include "CurrentData.h"
#include "Numeric.h"
#include "EigenPar.h"

/* Version commentee par A. Nicolet de Lanczos.c le 2001/11/29 */

/* Bibliographie

   Numerical Linear Algebra, Trefethen & Bau, SIAM, Philadelphia,
   1997.  Tres recent et tres clair, excellent ouvrage introductif sur
   l'algebre lineaire numerique moderne!
   
   Y. Saad, Numerical Methods for Large Eigenvalue Problems,
   Manchester University Press Avantages: gratuit sur le Web
   (http://www.users.cs.umn.edu/~saad/books.html et le bouquin sur la
   resolution des systemes s'y trouve aussi maintenant!), assez recent
   (1992), tres cible sur le probleme comme son titre l'indique, tres
   clair avec peu de prerequis (le premier chapitre vous rememore tout
   ce que vous devez savoir sur les matrices pour commencer), bon
   niveau theorique sans abstraction stratospherique et decrit
   precisement les ALGORITHMES (dont par exemple la methode
   d'Arnoldi-Tchebychev conu par Saad lui-meme).
   
   Golub & Van Loan, Matrix Computations, Johns Hopkins University
   Press, 3rd ed, 1996.  La troisieme edition de la bible du calcul
   matriciel applique!

   Valeurs propres des matrices, Franoise Chatelin, Masson, Paris,
   1988.  Bien cible et assez theorique, plutot court, un bon texte
   complementaire pour se faire un synthese.

   Analyse Numerique, sous la direction de Jacques Baranger, Hermann,
   Paris, 1991, Chapitre 7 par Franoise Chatelin. Bon chapitre
   synthetique par l'auteur de l'ouvrage precedent.

   M. Geradin, D. Rixen, Theorie des Vibrations, Masson, Paris, 2eme
   ed, 1996.  La reference utilisee par Benoit Meys. Introduction
   de la methode de Lanczos comme outil pour l'ingenieur donc approche
   tres pragmatique et tres instructive pour le passage aux
   applications.

   Isaacson & Keller, Analysis of Numerical Methods, Dover (edition
   originale 1966).  L'ouvrage general de reference en calcul
   numerique dans les annees 60.

   J. Stoer, R. Bulirsh, introduction to Numerical Analysis, Springer
   Verlag, 1979.  Un des meilleurs ouvrages generaux de reference en
   calcul numerique dans les annees 80. Entre Isaacson & Keller et
   Stoer & Bulirsh, il y a le Dahlquist et Bjork des annees 70 mais
   c'est un ouvrage beaucoup plus introductif.

   W.H. Press & al., Numerical Recipes, Cambridge University Press,
   ... Le best seller du calcul numerique des annees 90. Generation
   fast food oblige, c'est le fast programming. Indispensable en fait
   comme trousse de secours: un petit rappel theorique rapide ou une
   petite routine toute faite peuvent etre des gestes qui
   sauvent...

   Et enfin, les Wilkinson: J. H. Wilkinson, Rounding Errors in
   Algebraic Processes, Dover Petit ouvrage sur l'analyse d'erreur
   comme son nom l'indique.

   J. H. Wilkinson, C. Reinsch, Linear Algebra, Springer Verlag, 1971.
   Un pre-Numerical Recipes plus pousse sur l'algebre lineaire et avec
   tous les algorithmes en Algol (Je vous parle d'un temps que les
   moins de vingt ans ne peuvent pas connaitre...).

   J. H. Wilkinson, Algebraic eigenvalue problem, Clarendon Press,
   1965.  Le livre sacre du calcul des valeurs propres... que je n'ai
   jamais eu en main. Il va falloir que je mette la main dessus...
*/


/* Quelques modifications envisagees dans le futur (par ordre
   approximatif de difficulte et probablement par ordre chronologique
   de realisation future...) :

   1) Normer les vecteurs propres (valeur absolue max d'un element =
   1) -> OK mais ne resout pas entierement le probleme de l'affichage
   dans le post.

   2) Examiner de plus pres et eventuellement valider ou ameliorer
   quelques details comme par exemple je verrais bien le test if
   (fabs(shift) > 1.e-10) en relatif plutot qu'en absolu. De toute
   facon, si on essaye d'imposer un shift si petit, il vaudrait mieux
   afficher un message d'avertissement!  CHANGER LES NOMS DE CERTAINES
   VARIABLES

   3) Faire fonctionner LinAlg correctement pour les complexes (au
   moins le produit scalaire!!). En gros le codage des matrices et des
   vecteurs complexes du systeme est

   | Re Im |   et  | Re |
   |-Im Re |       | Im |

   Les matrices rectangulaires orthogonales sont stockees dans une
   liste de vecteurs Lan (Changer le nom en Q ou V ?).  Ce qui doit
   tre fait est le passage a une matrice de Hessenberg complexe Hes ->
   HesR,HesI ... mais on ne peut pas utiliser le codage ci-dessus car
   on perdrait la structure de Hessenberg.

   4) Calculer le RESIDU associe a un couple : la theorie permet de
   calculer facilement ce resultat, voir plus loin -> OK utilise pour
   selectionner les valeurs propres qui ont converge.

   5) Ameliorer les valeurs propres et les vecteurs propres a
   posteriori (il y a une methode simple expliquee dans Numerical
   Recipes, je crois...) et eventuellemnt introduire un critere
   d'arret.

   6) Amelioration de l'algorithme: Preconditionner (il y a des trucs
   la dessus dans Saad, je crois).  Reorthogonaliser les colonnes de
   Qn -> OK mais n'apporte pas grand chose !  Faire des redemarrages
   ... a essayer, ce n'est pas un gros travail, par exemple prendre
   comme point de depart une combinaison lineaire (somme !)  des
   vecteurs (normes !) qui ont deja plus ou moins converges ! Hn
   incomplet etc. (cf. Golub & Van Loan, Chatelin) L'algorithme de
   Saad 'Arnoldi-Tchebyshev' semble un bon candidat !

   7) 'Fixer omega et chercher gamma' conduit a un probleme aux VP
   generalise du type: lambda^2 M1 v + lambda M2 v + M3 v = 0 On peut
   mettre ce probleme sous la forme d'un probleme generalise de taille
   'double' du type

   | M1 M2 | |v1| lambda +  | 0  M3 | |v1| = 0
   | 0  I  | |v2|           |-I  0  | |v2|

   Comme d'habitude, on n'a besoin que du produit Matrice Vecteur ce
   qui revient a travailler avec les 'demi-vecteurs' v1, v2 de v et
   les sous-matrices M1, M2, M3: on doit calculer le vecteur 'double':

   | -M1 v2 + M2 M3^-1 v1 | 
   | M3^-1 v1             | 

   Il sera surtout necessaire d'avoir acces aux sous-matrices M1, M2,
   M3.  On peut egalement envisager le probleme des courbes de
   dispersion: Trouver des paires de nombres (omega,gamma) telles
   qu'il existe un vecteur v verifiant: gamma^2 M1 v + gamma M2 v +
   omega^2 M3 v + M4 v = 0 Le jeu consiste a faire varier (pas trop
   vite) omega et calculer les nouveaux gammas en utilisant si
   possible l'information fournie par les resolutions precedentes ...
   
   8) Autres methodes(puissance inverse?)
   
   Question: et si on utilisait une librairie comme ARPACK et son
   interface objet en C++ ARPACK++?
   http://www.caam.rice.edu/software/ARPACK/
*/

#if !defined(HAVE_GSL)

#include "nrutil.h"

#else

#include <gsl/gsl_rng.h>

/* base-1 index: numerical recipes convention */
#define dvector(dummy, A) (double *)Malloc(sizeof(double) * ((A) + 1))
#define dmatrix(dummy, A, dummy2, B) newmatrix((A) + 1, (B) + 1)
#define free_dvector(A, dummy, dummy2) Free(A)
#define free_dmatrix(M, dummy, A, dummy2, dummy3) freematrix((M), (A))

double **newmatrix(int nrow, int ncol){ /* accessible as a[i][j] */
  int i;
  double **a;
 
  a = (double **)Malloc(sizeof(double *) * nrow);
  for (i = 0; i < nrow; i++)
    a[i] = (double *)Malloc(sizeof(double) * ncol);
  return a;
}

void freematrix(double **a, int nrow){
  int i;
  double **b;
  
  b = a;
  for (i = 0; i < nrow; i++) Free(*b++);
  Free(a);
}

double gsl_random(){ /* uniform random numbers in range [0.0, 1.0) */
  double u;
  static const gsl_rng_type *T = NULL;
  static gsl_rng *r;
  
  if (T == NULL) { /* initialize */
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);
  }

  /* we never free this: gsl_rng_free(r); */

  u = gsl_rng_uniform(r);
  return u;
}

/* calc complex eigenvalues of Hessenberg form of real non-symmetrical
   matrix: GSL doesn't have it, so we use LAPACK's DHSEQR routine */

#if !defined(HAVE_BLAS_LAPACK)

void hqr(double **hin, int nhess, double *wrout, double *wiout){
  Msg(GERROR, "Lanczos algorithm not available without BLAS and LAPACK");
}

#else

#include "Compat.h"

#if defined(HAVE_UNDERSCORE)
#define dhseqr_ dhseqr
#endif

EXTERN_C_BEGIN
void dhseqr_(char *JOB, char *COMPZ, int *N, int *ILO, int *IHI, double *H, 
	     int *LDH, double *WR, double *WI, double **Z, int *LDZ,
	     double *WORK, int *LWORK, int *IRET);
EXTERN_C_END

void hqr(double **hin, int nhess, double *wrout, double *wiout){
  int n, ilo, ihi, ldh, ldz, lwork, info, i, j, k;
  double *h, *wr, *wi, *work, **dummy;
  char job[] = "E", compz[] = "N";
  
  n = nhess;
  ilo = 1;
  ihi = n;
  ldh = n;
  ldz = n;
  lwork = n;
  
  wr = (double *)Malloc(sizeof(double) * n);
  wi = (double *)Malloc(sizeof(double) * n);
  work = (double *)Malloc(sizeof(double) * n);
  
  h = (double *)Malloc(sizeof(double) * (n * n));
  
  /* C to F77 interface: Fortran uses "column-major ordering" in
     arrays, so transpose matrix and pack into h[] (or h[][]).  This
     may or may not work with other systems.  "cfortran.h" would be a
     better way */
  k = 0;
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
      h[k++] = hin[j + 1][i + 1]; /* hin: base-1 index */
    }
  }
  
  dhseqr_(job, compz, &n, &ilo, &ihi, h, &ldh, wr, wi, dummy, &ldz, work, 
	  &lwork, &info);
  if (info != 0)
    Msg(GERROR, "Lanczos/lapack dhseqr error = %d", info);
  
  for (i = 0; i < n; i++) {
    wrout[i + 1] = wr[i]; /* base-1 */
    wiout[i + 1] = wi[i];
  }
}