lanczos.c

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

#endif

#endif /* if !HAVE_GSL */

/* calcul des vecteurs propres d'une matrice de Hessenberg reelle
   (**T) de taille N dont on donne la valeur propre valp resultat en
   sortie dans *v */

void cal_vec_pr_T(double **T, int N, double valp, double *v){
  int      i,j,k;
  double **mat,norm=0.0;

  GetDP_Begin("cal_vec_pr_T");

  /* cette procedure devra etre reecrite pour une matrice de
     Hessenberg COMPLEXE */

  mat = dmatrix(1,N,1,N);
  
  for(i=1;i<=N;i++){
    for(j=1;j<=N;j++){
      if(i==j)
	mat[i][j]=T[i][j]-valp;
      else
	mat[i][j]=T[i][j];
    }
  }
  
  /* totalement instable si mat[][] est tres petit.  a
     changer. probleme a creuser d'une facon generale, considerer le
     raffinement des valeurs et vecteurs propres */

  /* resolution de (T - valp I) v = 0 systeme indetermine car v n'est
     defini qu'a une constante multiplicative pres on fixe donc v[N] a
     1 !!! pas toujours possible !!!!!  et on resout par substitution
     arriere grace a la forme Hessenberg */
   
  v[N]=1.0;
  for(i=N;i>1;i--){
    v[i-1]=0.0;
    for(j=N;j>=i;j--) v[i-1]-=mat[i][j]*v[j];
    if(mat[i][i-1] != 0.0)
      v[i-1]/=mat[i][i-1];
    else {
      Msg(INFO, " --- INVARIANT SUBSPACE FOUND ! --- ");
      /* verifier que la manoeuvre de sortie est valide !!! */
      for(k=i;k<=N;k++)	v[k]=0.0;
      v[i-1]=1.0;
    }
  }

  /* fixer la norme L2 de v a 1 */
  for(i=1;i<=N;i++) norm+=v[i]*v[i];
  norm = sqrt(norm);
  for(i=1;i<=N;i++) v[i]/=norm;

  GetDP_End ;
}

/* Algorithme de Lanczos (IL S'AGIT EN FAIT d' A R N O L D I POUR LES
   SYSTEMES NON SYMETRIQUES)

   sont transmis :
   - le pointeur d'un problemes avec ses matrices DofData_P
   - le nombre d'iterations a effectuer LanSize
   - une liste donnant les indexes des vecteurs propres a sauvegarder
   - le decallage shift
*/

void Lanczos (struct DofData * DofData_P, int LanSize, List_T *LanSave, double shift){
  gVector  *Lan, *b, *x ;
  gMatrix  *K, *M ;
  int      i, j, k, ii, jj, NbrDof, Restart, ivmax, newsol ; 
  double   dum, dum1, dum2;
  double   **Hes, **IMatrix, *diag, *wr, *wi, *err ;
#if !defined(HAVE_GSL)
  long     mun = -1 ;
#endif
  struct Solution Solution_S ;
  struct EigenPar eigenpar;

  GetDP_Begin("Lanczos");

  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

  Msg(BIGINFO, "Lanczos - December 2001 - beta 0.2 A. Nicolet - Marseille ");

  if(!DofData_P->Flag_Init[1] || !DofData_P->Flag_Init[3])
    Msg(GERROR, "No System available for Lanczos: check 'DtDt' and 'GenerateSeparate'") ;

  /* lecture des parametres dans le fichier 'eigen.par' */
  EigenPar("eigen.par", &eigenpar);
  eigenpar.size = LanSize; /* Hack: we force this */

  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

  NbrDof = DofData_P->NbrDof ; /* taille du systeme */

  /* reservation de LanSize+1 vecteurs de taille NbrDof reperes par
     les adresses &Lan[i] dans la suite &Lan[i] sera note q_i dans les
     commentaires (colonnes d'une matrice rectangulaire
     orthogonale) */
  Lan = (gVector*)Malloc((LanSize+1)*sizeof(gVector)) ;
  for (i=0 ; i<LanSize+1 ; i++) 
    LinAlg_CreateVector(&Lan[i], &DofData_P->Solver, NbrDof,
			DofData_P->NbrPart, DofData_P->Part);

  /* resolution du PVP generalise K v = valp M v */
  /* identifier les matrices avec des matrices du probleme en cours */
  K = &DofData_P->M1 ; /* matrice des autres termes */
  /* L = &DofData_P->M2 ; */ /* matrice des termes en Dt pour le futur */
  M = &DofData_P->M3 ; /* matrice des termes en DtDt  */
  b = &DofData_P->b  ; 
  x = &DofData_P->CurrentSolution->x ;

  /* Si on veut traiter le probleme 'omega fixe' en propagation, les
     valeurs propres sont dans ce cas associees a la constante de
     propagation et donc a la derivation en z! La notation Dt est donc
     mal choisie dans ce cas ... */
  
  /* Cf. egalement remarque sur les courbes de dispersion! */
  
  /* On est ici dans un cas tres particulier ou tout est complexe sauf
     la matrice de Hessenberg et les valeurs propres. On fait Arnoldi
     car on construit une matrice de Hessenberg mais c'est en fait une
     matrice tridiagonale (-> beaucoup de calcul pour rien !  mais a
     ne ralenti pas trop le calcul car ici c'est le NOMBRE DE
     RESOLUTIONS DU SYSTEME pour les iterations inverses qui est
     couteux). Par contre on ne considere que la partie reelle des
     produits scalaires de vecteurs et on construit une matrice de
     Hessenberg reelle -> Ce n'est pas assez general pour etre du
     vrai Arnoldi sur une matrice quelconque !!!  */
  
  /* Construire une Hessenberg complexe !! */
  
  /* declaration des espaces de travail */
  diag    = dvector(1, LanSize+1);
  wr      = dvector(1, LanSize+1);
  wi      = dvector(1, LanSize+1);
  err     = dvector(1, LanSize+1);
  IMatrix = dmatrix(1, LanSize+1, 1, LanSize+1);
  Hes     = dmatrix(1, LanSize+1, 1, LanSize+1);

  /* initial random vector b=q_o pas optimal pour la reproductibilite
     des resultats ! ! !  pourquoi ne pas essayer des 1 partout ?
     Arnoldi-Tchebychev consiste a ameliorer le vecteur de depart
     d'Arnoldi a l'aide de Tchebychev */
#if defined(HAVE_GSL)
  for (i = 0; i < NbrDof; i++)
    LinAlg_SetDoubleInVector(gsl_random(), &Lan[0], i);
#else
  for (i=0 ; i<NbrDof ; i++) 
    LinAlg_SetDoubleInVector(ran3(&mun), &Lan[0], i) ;
#endif

  /* shifting: K = K - shift * M */
  /* DANGER DANGER depend de la mise a l'echelle de K et M */
  if (fabs(shift) > 1.e-10)
    LinAlg_AddMatrixProdMatrixDouble(K, M, -shift, K) ; 

  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

  /* iterations d' A R N O L D I */

  /* Soit le probleme au vp : A v = lambda v avec A matrice m * m et
     soit un vecteur arbitraire b (choix ????)  Kn la matrice m*n de
     Krylov donnee par <b, Ab, A^2 b, ... A^(n-1)b> (espace de Krylov
     K_n(A,b) = espace genere par la combinaison lineaire de ces
     vecteurs = lin(b, Ab, A^2 b, ... A^(n-1)b)) A l'etape n des
     iterations d'Arnoldi, la matrice m*n orthogonale Qn est le
     facteur de la decomposition QR de Kn = Qn Rn (m*n = (m*n) *
     (n*n)) La matrice de Hessenberg Hn correspond a la projection Hn
     = Qn^* A Qn (n*n = (n*m) * (m*m) * (m*n) ) (^* conjugue hermitien
     ou adjoint - l'algorithme est valable pour les matrices COMPLEXES
     NON HERMITIENNES ) Les iterations successives sont reliees par A
     Qn = Qn+1 H'n ou H'n est la matrice (n+1)*n obtenue en completant
     Hn avec l'element h_(n+1,n).

     La structure des iterations est (A est une matrice, b,v des
     vecteurs colonnes et les q des vecteurs colonnes de Qn, les h
     sont les coefficients de Hn)

     b arbitraire, q1 = b/norme(b)
     POUR n=1,2,3,...
          |v = A q_n
          |POUR j = 1,n 
          |     |h_(j,n) = q_j^* v
          |     |v = v - h_(j,n) q_j
          |h_(n+1,n)=norme(v)               
          |q_(n+1) = v/h_(n+1,n)     

    (Que se passe-t'il si h_(n+1,n)=0 ? -> on a trouve un sous espace
    invariant = un 'sous-espace propre')

    Les valeurs propres sont approximees par les valeurs propres de Hn
    dans le sens suivant : Une matrice A verifie son polynome
    caracteristique P(z) soit P(A) = 0 soit norme(P(A) b) = 0 pour
    tout vecteur b Soit pn(z) le polynome caracteristique de Hn
    (polynome d'Arnoldi) C'est le polynome de degre n qui minimise
    norme(p(A) b) pour tout polynome p(z) de degre n.  Remarque : le
    polynome d'Arnoldi ideal ou polynome de Tchebychev de A est le
    polynome p_Tcheb(z) de degre n qui minimise norme(p(A)) (norme
    matricielle) Ce polynome ne depend pas d'un vecteur b arbitraire
    mais il n'existe pas d'algorithme efficace pour le calculer

    Au cours des iterations on construit Hn et Qn : Hn = Qn* A Qn avec
    Qn* Qn = I et on reecrit A Qn = Qn+1 H'n comme A Qn = Qn Hn +
    h_(n+1,n) q_n+1 e^*_n (h_n+1,n prochain coefficient de Hn+1, q_n+1
    prochaine colonne de Qn+1, e_n base canonique de C^n)

    Si v est vecteur propre de Hn alors Hn v = lambda v.  On choisit
    Norme(v)=1.  v est un vecteur de Ritz et lambda une valeur de
    Ritz.  Ils constituent l'approximation au sens de Galerkine dans
    l'espace de Krylov K_n(A,b) du probleme aux valeurs propres : <w,
    A v - lambda v> = 0 pour tout w dans K_n(A,b)

    Les coefficients de h sont les coefficients du developpement de A
    dans la base orthogonale de K_n(A,b) constituees des colonnes de
    Qn
     
    Donc Qn* A Qn v = lambda v Malheureusement Qn n'est pas carree et
    ce n'est pas une relation de similitude.  Neanmoins si on poursuit
    le processus jusque n=m (si possible) Qm diagonalise A et on a une
    relation de similitude.  Dans le cas d'une relation de similitude
    si v est vecteur propre de Q* A Q alors Q* A Q v = lambda v et
    comme Q Q* = I on a A Q v = lambda Q v et Q v est donc vecteur
    propre de A.  Ici on a Qn* Qn = I(n*n) mais pas Qn Qn* = I(m*m) Qn
    Qn* est le projecteur orthogonal de l'espace C^m (dimension m) (=
    vecteurs quelconques de m complexes) sur l'espace de Krylov
    K_n(A,b) (dimension n).  On estime le vecteur propre de A par un
    'relevement' de v a l'aide de Qn: (estimation du vecteur propre de
    A) = Qn v.  La multiplication par Qn est donc un peu l'equivalent
    de l'interpolation.

    On obtient une approximation du residu de la facon suivante On
    pose x = Qn v Norme(A x - lambda x)= Norme(A Qn v - lambda Qn v) =
    Norme(A Qn v - Qn Hn v). Norme((A Qn - Qn Hn) v) = Norme(h_(n+1,n)
    q_n+1 e^*_n v) =h_(n+1,n) v(n) car Norme(q_n+1)=1 =dernier coef de
    Hessenberg x derniere composante du vecteur de Ritz qu'il suffit
    de comparer a lambda pour avoir une estimation de l'erreur
    relative !

    Probleme generalise : K v = valp M v -> utiliser le produit
    scalaire <M x, y> ou M joue le role de metrique.  Arnoldi/Lanczos
    donne les plus grandes valeurs propres -> utiliser la matrice
    inverse K*-1 -> on aboutit a une resolution de systeme.

    A venir, les PVP du second ordre du type (M1 valp^2 + M2 valp +
    M3) v = 0

    Il faudra egalement voir si on ne peut pas ameliorer le processus
    en utilisant un PRECONDITIONNEMENT pour les problemes aux VP
    (different de ceux pour la resolution des systemes) Voir Y. Saad,
    Numerical Methods for Large Eigenvalue Problems disponible sur le
    Web un redemarrage implicite, voir doc ARPACK
      
    NOTATION du PROGRAMME
    Dans ce qui suit K et M portent ce nom
    Hes contient Hn ( Hes[i][j]=Hn_(i,j) )
    &Lan[i] est l'adresse de la colonne q_i
    Variable de boucle exterieure i a la place de n
                       interieure j a la place de i
    L'indice i commence a zero. */


  /* traitement special des premieres iterations */
  /* ------------------------------------------- */

  /* etape 0 */

  /* b = M q_o */
  LinAlg_ProdMatrixVector(M, &Lan[0], b);
  /* dum = b^T q_o  (^T = transposition) */
  LinAlg_ProdVectorVector(b, &Lan[0], &dum);
  /* Lan[0] is built: q_o = q_o/sqrt(dum) */
  LinAlg_ProdVectorDouble(&Lan[0], 1./sqrt(dum), &Lan[0]);

  Msg(INFO, "Lanczos iteration 1/%d", LanSize);

  /* etape 1 */

  /* b = M * Lan[0]: b = M q_o */
  LinAlg_ProdMatrixVector(M, &Lan[0], b); 
  /* Lan[1] = K^-1 * b: q_1 = K^-1 b */  
  LinAlg_Solve(K, b, &DofData_P->Solver, &Lan[1]);

  /* pas d'inversion explicite, on utilise le solver bien sur ! -> on
     a calcule K^-1 M q_o */