lanczos.c

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

  /* alpha1 = Lan[0]^T * M * Lan[1]: dum1 = b^T q_1*/
  LinAlg_ProdVectorVector(b, &Lan[1], &dum1); 

  /* orthogonalization: Lan[1] -= alpha1 * Lan[0]: q_1 = q_1 - dum1 q_o */
  LinAlg_AddVectorProdVectorDouble(&Lan[1], &Lan[0], -dum1, &Lan[1]); 

  /* b = M * Lan[1] in prevision: b = M q_1 */
  LinAlg_ProdMatrixVector(M, &Lan[1], b); 
  /* gama2 = beta1 = sqrt(Lan[1]^T * M * Lan[1]): dum2 = b^T q_1 */
  LinAlg_ProdVectorVector(b, &Lan[1], &dum2);

  /* normation in the metric M: Lan[1] = q_1  is built */  
  dum2 = sqrt(dum2); 
  LinAlg_ProdVectorDouble(&Lan[1], 1./dum2, &Lan[1]); 
  /* b = M * Lan[1] in prevision of next step */
  LinAlg_ProdVectorDouble(b, 1./dum2, b);

  /* debut de construction de la matrice de Hessenberg pour
     l'approximation des vp */
  Hes[1][1] = dum1;
  Hes[2][1] = dum2;

  /* print_lanczos_to_file (1, NbrDof, Hes, Lan, shift, Name); */

  /* Debut de la boucle principale des iterations d' A R N O L D I */
  /* ------------------------------------------------------------- */

  Restart = 2 ; /* pour 'Arnoldi with restarting', cf. Golub & Van Loan */
  for (i=Restart ; i<=LanSize ; i++){
    
    Msg(INFO, "Lanczos iteration %d/%d", i, LanSize);

    /* q_1 = K^-1 b  avec b deja multiplie par M */
    LinAlg_Solve(K, b, &DofData_P->Solver, &Lan[i]);
    
    for (j=1 ; j<=i ; j++){
      /* x = M q_(j-1) */
      LinAlg_ProdMatrixVector(M, &Lan[j-1], x);
      /* h_(j,i) = x^T q_i */
      LinAlg_ProdVectorVector(x, &Lan[i], &Hes[j][i]);
      /* orthogonalization: q_i = q_i - h_(j,i) q_(j-1) */
      LinAlg_AddVectorProdVectorDouble(&Lan[i], &Lan[j-1], -Hes[j][i], &Lan[i]);
    }

    /* OPTIONNAL PART : re-orthogonalization DGKS */
    if (eigenpar.reortho == 1) {
      Msg(INFO, " *** reorthogonalization *** ");
      for (j=1 ; j<=i ; j++) {
	LinAlg_ProdMatrixVector(M, &Lan[j-1], x);
	LinAlg_ProdVectorVector(x, &Lan[i], &dum);
	diag[j]=dum;
      }
      /* two separate loops so that &Lan[i] is not modified in the first one */
      for (j=1 ; j<=i ; j++) {
	/* reorthogonalization */ 
	LinAlg_AddVectorProdVectorDouble(&Lan[i], &Lan[j-1], -diag[j], &Lan[i]);
	Hes[j][i]=Hes[j][i]+diag[j];
	Msg(INFO, " hes %d = %g corrected by %g ",j,Hes[j][i],diag[j]);
      }
    }

    /* x = M q_i */
    LinAlg_ProdMatrixVector(M, &Lan[i], x);
    /* dum  = x^T q_i */
    LinAlg_ProdVectorVector(x, &Lan[i], &dum);
    Hes[i+1][i] = sqrt(dum); 
    /* q_i = q_i / h_(j,i)*/ 
    LinAlg_ProdVectorDouble(&Lan[i], 1./Hes[i+1][i], &Lan[i]);
    /* b = M * Lan[i] in prevision of next step */
    LinAlg_ProdMatrixVector(M, &Lan[i], b);
    
    /* eigen value computation of the current Hes matrix (On complete
       la matrice de Hessenberg par des zeros) */
    for(ii=3 ; ii<=i ; ii++)
      for(jj=1 ; jj<=ii-2 ; jj++)
	Hes[ii][jj] = 0.0 ;
    
    for(ii=1 ; ii<=i ; ii++)
      for(jj=1 ; jj<=i ; jj++)
	IMatrix[ii][jj] = Hes[ii][jj] ;

    /* cette partie est pour info, elle n'est necessaire que pour
       afficher les VP on pourrait imaginer de ne pas la calculer
       chaque fois! */
    hqr(IMatrix, i, wr, wi) ;

    /* Algorithme QR pour une matrice de Hessenberg (from Numerical
       Recipes) Pour une MATRICE REELLE ???????????????????  Ecrire
       l'algorithme correspondant en complexe : decomposition LR cf
       version Algol dans Wilkinson */

#if 0 
    print_valpr_to_file(i,wr,wi,shift,Name);
    print_lanczos_to_file (i,NbrDof,Hes,Lan,shift,Name);

    Msg(INFO, "Lanczos eigenvalue of the transformed problem ");
    for(k=1 ; k<=i ; k++)
      Msg(INFO, "Lanczos eigenvalue %d = %g +I %g",k, 
	  wr[k], wi[k]);

    Msg(INFO, "Lanczos eigenvalue of the original problem ");
    for(k=1 ; k<=i ; k++)
      Msg(INFO, "Lanczos eigenvalue %d = %g (%g) %g",k, 
	  shift+1.0/wr[k], sqrt(shift+1.0/wr[k])/TWO_PI, wi[k]);
#endif

    /* ATTENTION shift+1.0/wr[k] ne donne la VP reelle que si wr[i]
       est NEGLIGEABLE. Quid des VP complexes ? */

    /* store the real eigen values */
    for (k=1 ; k<=i ; k++)
      wi[k] = shift+1.0/wr[k];

    /* estimation d'erreur et test de convergence */
    Msg(INFO, "------------------ hessenberg coeff %d = %g  ",i,Hes[i+1][i]);
    if (Hes[i+1][i] < 1e-20) 
      Msg(INFO, " --- INVARIANT SUBSPACE FOUND ! --- ");
    
    /* search the largest eigenvalue */
    dum = 0. ; ivmax =1 ;
    for (k=1 ; k<=i ; k++)
      if (wr[k] > dum) {ivmax = k; dum = wr[k];};
    
    if (wr[ivmax] == 0.) Msg(WARNING," OOOPS !! - no positive eigen value ? - ");
    
    Msg(INFO, "Max eigenvalue = %g on %d ", dum, ivmax);  

    /* compute the corresponding eigenvector */
    cal_vec_pr_T(Hes, i, wr[ivmax], diag);
    Msg(INFO, "Last eigenvector component = %g ", diag[i]);  
    Msg(INFO, " ******** Residual estimate = %g ",
	Hes[i+1][i] * diag[i] / wr[ivmax] );
  }

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

  Msg(INFO, "Final eigenvalue/eigenvector Computation");

  /* eigen value computation of the final Hes matrix (Calcul final des
     valeurs propres = Euh ? Deja fait en sortie de boucle principale,
     non ??) */

  for(ii=3 ; ii<=LanSize ; ii++)
    for(jj=1 ; jj<=ii-2 ; jj++)
      Hes[ii][jj] = 0.0 ;

  for(ii=1 ; ii<=LanSize ; ii++)
    for(jj=1 ; jj<=LanSize ; jj++)
      IMatrix[ii][jj] = Hes[ii][jj] ;

  hqr(IMatrix, LanSize, wr, wi);

  /* store the real eigen values (est-ce une bonne idee de stocker les
     VP reelles reconstituees dans la partie complexe ???? non!) */
  for (k=1 ; k<=LanSize ; k++)
    wi[k] = shift+1./wr[k]; 

  /* eigen vector computation of the final Hes matrix (Pour chacune
     des valeurs propres de Hn, calcul du vecteur propre a l'aide de
     la routine 'cal_vec_pr_T') */
  for(i=1 ; i<=LanSize ; i++){
    /* diag ne sert qu'ici, son nom est loufoque ! */
    cal_vec_pr_T(Hes, LanSize, wr[i], diag);
    
    /* estimation d'erreur et test de convergence */
    /* Msg(INFO, "*********** estim %d = %g", i, 
       Hes[LanSize+1][LanSize]*diag[LanSize]); */
    if (wr[i]>1e-20) 
      err[i] = Hes[LanSize+1][LanSize]*diag[LanSize]/wr[i];
    else
      err[i] = 1.e99; /* a changer! */

    /* copy the current eigen vector in IMatrix */
    for(j=1 ; j<=LanSize ; j++)
      IMatrix[j][i]=diag[j];
  }
  
  /* reconstruction of the global eigen vectors */
  
  newsol = 0;

  for(i=0 ; i<List_Nbr(LanSave) ; i++){
    List_Read(LanSave, i, &ii) ;
    
    if(ii<1 || ii>LanSize){
      Msg(WARNING, "Eigenvalue index out of range") ;
      break ;
    }
    
    /* test de validite de la valeur propre */
    Msg(BIGINFO, "Eigenvalue %d = %.16g (f = %.16g)",
	ii, wi[ii], sqrt(wi[ii])/TWO_PI);

    /* if error smaller than required precision and non pathological
       eigenvalue ! */
    if ((err[ii] < eigenpar.prec ) && (wr[ii]>0.)) { 
      Msg(BIGINFO, "GOOD -> Eigenvalue %d = %g with error %g ",
	  ii, wr[ii], err[ii]);
      
      if (newsol) {
	/* create new solution */
	LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, NbrDof,
			    DofData_P->NbrPart, DofData_P->Part);
	List_Add(DofData_P->Solutions, &Solution_S) ;
	DofData_P->CurrentSolution = (struct Solution*)
	  List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1) ;
      } 
      newsol = 1;

      /* Time = partie reelle de la vp ? recuperation possible dans le
         Post ? J'ai verifie dans l'affichage de Gmsh et a semble bien
         tre le cas */
      DofData_P->CurrentSolution->Time = wi[ii] ;
      DofData_P->CurrentSolution->TimeImag = 0. ;
      DofData_P->CurrentSolution->TimeStep = (int)Current.TimeStep ;
      DofData_P->CurrentSolution->TimeFunctionValues = NULL ;
      DofData_P->CurrentSolution->SolutionExist = 1 ;
      LinAlg_ZeroVector(&DofData_P->CurrentSolution->x) ;

      /* increment the global timestep counter so that a future
	 GenerateSystem knows which solutions exist */
      Current.TimeStep += 1.;
      /* update the current value of Time and TimeImag so that
	 $EigenvalueReal and $EigenvalueImag are up-to-date */
      Current.Time = wi[ii];
      Current.TimeImag = 0.;

      /* boucle de taille m = taille des matrices A,K,M du probleme */
      /* calcul de la composant k du vecteur ii */
      for(k=0; k<NbrDof; k++){
	
	/* boucle de taille n = le nombre de vecteurs propres a estimer */
	/* produit vecteur^T vecteur mais PAS avec un vecteur q */
	/* on balaye Qn dans l'autre sens ! */
	for (j=1 ; j<=LanSize ; j++){
	  /* dum = element k de q_(j-1) */	  
	  LinAlg_GetDoubleInVector(&dum, &Lan[j-1], k) ;

	  /* x_k = x_k + q_(k,j-1)* v_j */
	  /* on a multiplie le 'ii ieme' vecteur propre de Hn par Qn */  
	  /* v de taille n -> Qn v de taille (m*n) * n = m     */
	  LinAlg_AddDoubleInVector(IMatrix[j][ii]*dum, 
				   &DofData_P->CurrentSolution->x, k) ;
	}
      }
      
      /* normation L infini : abs plus grand element mis a un */
      /* Msg(INFO, "normation du vecteur propre %d  ",ii); */
      /* determination du maximum */
      dum = 0.; dum1 = 0.;
      for(k=0; k<NbrDof; k++){
	LinAlg_GetDoubleInVector(&dum,&DofData_P->CurrentSolution->x, k);
	if (fabs(dum)>dum1)
	  dum1 = dum;
      }
      /* division par le max */
      if (dum1 > 1.e-16) 
	LinAlg_ProdVectorDouble(&DofData_P->CurrentSolution->x,1./dum1,
				&DofData_P->CurrentSolution->x);
      
      /* estimation d'erreur et test de convergence */
      /* Msg(INFO, "------------------ estim %d = %g  ",ii,Hes[ii+1][ii]); */
      
    } /* fin du test de validite */
    else{
      Msg(BIGINFO,"BAD! -> Eigenvalue %d = %g with error %g ", ii, wr[ii], err[ii]); 
    }
    
  }    
  
  /* c'est fini */
  
  /* tester newsol pour voir si au moins une solution a ete trouvee ! */
  
  /* desallocation */
  
  for (i=0 ; i<LanSize ; i++) LinAlg_DestroyVector(&Lan[i]);
  Free(Lan) ;
  
  free_dvector(diag, 1, LanSize);
  free_dvector(wr, 1, LanSize);
  free_dvector(wi, 1, LanSize);
  free_dvector(err, 1, LanSize);
  free_dmatrix(IMatrix, 1, LanSize, 1, LanSize);
  free_dmatrix(Hes, 1, LanSize, 1, LanSize);

  GetDP_End ;
}