matricecreuse_tpl.hpp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

  case MatriceElementaire<R>::Symmetric : //throwassert(L ==U);   
    for (il=0; il<me.n; ++il) // modif overflow FH win32
      { i=mi[il];
      for (jl=0;jl<= il;++jl)
       { j=mj[jl]  ;
	 if      (j<i)  L[ pL[i+1] - (i-j) ] += *al++;
	 else if (j>i)  U[ pU[j+1] - (j-i) ] += *al++;
	 else           D[i] += *al++;}}
    break;
  default:
    cerr << "Big bug type MatriceElementaire unknown" << (int) me.mtype << endl;
    exit(1);
    break; 
  }      
  return *this;
} 

template<class R>
ostream& MatriceProfile<R>::dump (ostream& f) const 
{f<< " matrix skyline " << this->n << '\t' << this->m << '\t' ;
 f <<  "  this " << endl;
 f << " pL = " << pL << " L ="  << L << endl
   << " pU = " << pU << " U ="  << U << endl
   << " D = " << D << endl;
 if ( (pL == pU) &&  (U == L) )
   if (pL && L) 
     {f << " skyline symmetric " <<endl;
     int i,j,k;
     for (i = 0;i<this->n;i++) 
       { f << i << " {" << pL[i+1]-pL[i] << "}" <<'\t' ;
       for (k=pL[i];k<pL[i+1];k++)
	 { j=i-(pL[i+1]-k);
	 f << j << " " << L[k] << "; "; 
	 }
       f <<  i  << ":" << D[i] << endl  ;
       }
     }
   else f << "Skyline: pointeur null " <<endl; 
 else 
   { 
     f << " Skyline  non symmetric " << endl;
     int i,k;
     for (i = 0;i<this->n;i++) 
       {
	 f << i ;
	 if (pL && L) 
	   {
	     f << " jO=" << i-pL[i+1]+pL[i] << " L= " <<'\t' ;
	     for (k=pL[i];k<pL[i+1];k++)
	       { 
		 f <<  " " << L[k] ; 
	       }
	   }
	 if (D)
	   f  << " D= " << D[i]  << '\t' ;
	 else 
	   f  << " D=0 => 1 ; ";
	 if (pU && U) 
	   {
	     f << " i0=" << i-pU[i+1]+pU[i] << " U= " <<'\t' ;
	     for (k=pU[i];k<pU[i+1];k++)
	       f << " " << U[k] ; 

	   }
	 f << endl;  
       }
    
   }
 return f;
}
template<class R>
void MatriceProfile<R>::cholesky(double eps) const {
  double eps2=eps*eps;
  R  *ij , *ii  , *ik , *jk , xii;
  int i,j,k;
  if (L != U) ERREUR(factorise,"Skyline matrix non symmetric");
  U = 0; // 
  typefac = FactorizationCholeski;
  if ( norm(D[0]) <= 1.0e-60)
      ERREUR(cholesky,"pivot (" << 0 << ")= " << D[0] )
  
  D[0] = sqrt(D[0]); 
  ij = L ; // pointeur sur le terme ij de la matrice avec j<i 
  for (i=1;i<this->n;i++) // boucle sur les lignes 
    { ii = L+pL[i+1]; // pointeur sur le terme fin de la ligne +1 =>  ij < ii;
    xii = D[i] ; 
    for ( ; ij < ii ; ij++) // pour les j la ligne i
      { j = i -(ii - ij); 
      k = Max( j - (pL[j+1]-pL[j]) ,  i-(pL[i+1]-pL[i]) ); 
      ik =  ii - (i - k); 
      jk =  L + pL[j+1] -(j - k); 
      k = j - k ;
      R s= -*ij; 
#ifdef WITHBLAS
      s += blas_sdot(k,ik,1,jk,1);
#else
      while(k--) s += *ik++ * *jk++;  
#endif
      *ij =  -s/D[j] ;
      xii -= *ij * *ij ;
      }
    // cout << norm(xii) << " " << Max(eps2*norm(D[i]),1.0e-60) << " " << sqrt(xii) <<endl;
    if ( norm(xii) <= Max(eps2*norm(D[i]),1.0e-60)) 
      ERREUR(cholesky,"pivot (" << i << ")= " << xii << " < " << eps*abs(D[i]))
    D[i] = sqrt(xii);
    }
}
template<class R>
void MatriceProfile<R>::crout(double eps) const  {
  R  *ij , *ii  , *ik , *jk , xii, *dkk;
  int i,j,k;
  double eps2=eps*eps;
  if (L != U) ERREUR(factorise,"Skyline matrix  non symmetric");
  U = 0; // 
  typefac = FactorizationCrout;
   
  ij = L ; // pointeur sur le terme ij de la matrice avec j<i 
  for (i=1;i<this->n;i++) // boucle sur les lignes 
    { ii = L+pL[i+1]; // pointeur sur le terme fin de la ligne +1 =>  ij < ii;
    xii = D[i] ; 
    for ( ; ij < ii ; ij++) // pour les j la ligne i
      { j = i -(ii - ij); 
      k = Max( j - (pL[j+1]-pL[j]) ,  i-(pL[i+1]-pL[i]) ); 
      ik =  ii - (i - k); 
      jk =  L + pL[j+1] -(j - k); 
      dkk = D + k;
      k = j - k ; 
      R s=-*ij;
      while ( k-- ) s += *ik++ * *jk++ * *dkk++;  
      *ij = -s/ *dkk ; // k = j ici 

      xii -= *ij * *ij * *dkk;
      }
    if (norm(xii) <= Max(eps2*norm(D[i]),1.0e-60))
      ERREUR(crout,"pivot (" << i << " )= " << abs(xii)<< " <= " << eps*abs(D[i]) << " eps = " << eps)
	D[i] = xii;
    }
}
template<class R>
void MatriceProfile<R>::LU(double eps) const  {
  R s,uii;
  double eps2=eps*eps;
  int i,j;
  if (L == U && ( pL[this->n]  || pU[this->n] ) ) ERREUR(LU,"matrix LU  symmetric");
  if(verbosity>3)
  cout << " -- LU " << endl;
  typefac=FactorizationLU;

  for (i=1;i<this->n;i++) // boucle sur les sous matrice de rang i 
    { 
      // for L(i,j)  j=j0,i-1
      int j0 = i-(pL[i+1]-pL[i]);
      for ( j = j0; j<i;j++)
        {           
          int k0 = Max(j0,j-(pU[j+1]-pU[j]));
          R *Lik = L + pL[i+1]-i+k0; // lower
          R *Ukj = U + pU[j+1]-j+k0; // upper
          s =0;
#ifdef WITHBLAS
	  s = blas_sdot(j-k0,Lik,1,Ukj,1);
	  Lik += j-k0;
#else
          for (int k=k0;k<j;k++) // k < j < i ;
	    s += *Lik++ * *Ukj++ ;     // a(i,k)*a(k,j);
#endif
	  *Lik -= s;
	  *Lik /= D[j]; //  k == j here
        }
      // for U(j,i) j=0,i-1        
      j0=i-pU[i+1]+pU[i];
      for (j=j0;j<i;j++) 
	{
	  s = 0;
	  int k0 = Max(j0,j-pL[j+1]+pL[j]);
	  R *Ljk = L + pL[j+1]-j+k0;   
	  R *Uki = U + pU[i+1]-i+k0;   
#ifdef WITHBLAS
	  s = blas_sdot(j-k0,Ljk,1,Uki,1);
          Uki += j-k0;
#else
	  for (int k=k0  ;k<j;k++)    // 
	    s +=  *Ljk++ * *Uki++ ;
#endif
	  *Uki -= s;  // k = j here 
	}
      // for D (i,i) in last because we need L(i,k) and U(k,i) for k<j
      int k0 = i-Min(pL[i+1]-pL[i],pU[i+1]-pU[i]);
      R *Lik = L + pL[i+1]-i+k0; // lower
      R *Uki = U + pU[i+1]-i+k0; // upper
      s =0;
#ifdef WITHBLAS
       s = blas_sdot(i-k0,Lik,1,Uki,1);
#else
      for (int k=k0;k<i;k++) // k < i < i ;
	s += *Lik++ * *Uki++ ;     // a(i,k)*a(k,i);
#endif
      // cout << " k0 " << k0 << " i = " << i << " " <<  s << endl;
      uii = D[i] -s;
      
      if (norm(uii) <= Max(eps2*norm(D[i]),1.0e-30))
	ERREUR(LU,"pivot (" << i << " )= " << abs(uii) << " <= " << eps*abs(D[i]) << " eps = " << eps);     
      
      D[i] = uii;
      
    }
}


template<class R>
KN_<R> & operator/=(KN_<R> & x ,const MatriceProfile<R> & a) 
{
  // --------------------------------------------------------------------
  //   si La diagonal D n'existe pas alors on suppose 1 dessus (cf crout)
  // --------------------------------------------------------------------
  R * v = &x[0];
  int n = a.n;  
  if (x.n != n ) 
    ERREUR (KN_<R> operator/(MatriceProfile<R>),"  matrice et KN_<R> incompatible");
  const R *ij ,*ii, *ik, *ki;
  R *xk,*xi;
  int i;
  switch (a.typefac) {
  case FactorizationNO:
    if (a.U && a.L) {cerr << "APROGRAMMER (KN_<R><R>::operator/MatriceProfile)";throw(ErrorExec("exit",2));}
   
    if ( a.U && !a.L ) 
      { // matrice triangulaire superieure
	// cout << " remonter " << (a.D ? "DU" : "U") << endl;
	ki = a.U + a.pU[n]; 
	i = n;
	while ( i-- )
	  { ii = a.U + a.pU[i];
          xi= xk  = v +  i ;
          if (a.D) *xi /= a.D[i];// pour crout ou LU
          while ( ki > ii) 
	    *--xk  -=  *--ki *  *xi ; 
	  }
      }
    else if  ( !a.U && a.L ) 
      { // matrice triangulaire inferieure
	// cout << " descente "  <<( a.D ? "LD" : "L" ) <<endl;
	ii = a.L;
	for (i=0; i<n; i++)
	  { ij = ik = (a.L + a.pL[i+1]) ;  // ii =debut,ij=fin+1 de la ligne 
          xk = v + i;
          R ss = v[i]; 
          while ( ik > ii) 
	    ss -= *--ik * *--xk ; 
          if ( a.D) ss /= a.D[i];// pour crout ou LU
          v[i] = ss ;
          ii = ij;
	  }
      }
    else if (a.D) 
      { // matrice diagonale
	// cout << " diagonal D" <<endl;
	for (i=0;i<n;i++) 
	  v[i]=v[i]/a.D[i];
      }
    break;
  case FactorizationCholeski:  
    //     cout << " FactorizationChosleski" << endl;
    x /= a.ld();
    x /= a.ldt();   
    break;
  case FactorizationCrout:
    //   cout << " FactorizationCrout" << endl;
    x /= a.l();
    x /= a.d();
    x /= a.lt();
    break;
  case FactorizationLU:
    //  cout << " FactorizationLU" << endl;
    x  /= a.l();
    x  /= a.du();
    break;
    /*   default:
	 ERREUR  (operator /=(MatriceProfile," Error unkown type of Factorization  =" << typefac);
    */
  }
  return x;
}

template <class R> 
 MatriceMorse<R>::MatriceMorse(KNM_<R> & A,double tol)
    :MatriceCreuse<R>(A.N(),A.M(),false),solver(0) 
      {
  double tol2=tol*tol;    
  symetrique = false;
  this->dummy=false;
  int nbcoeff=0;
  for(int i=0;i<this->n;i++)
    for(int j=0;j<this->m;j++)
      if(norm(A(i,j))>tol2) nbcoeff++;

  nbcoef=nbcoeff;
  nbcoeff=Max(nbcoeff,1); // pour toujours alloue quelque chose FH Bug dans CheckPtr
  a=new R[nbcoeff] ;
  lg=new int [this->n+1];
  cl=new int [nbcoeff];
  nbcoeff=0;
  R aij;
  for(int i=0;i<this->n;i++)
   { 
    lg[i]=nbcoeff;
    for(int j=0;j<this->m;j++)
     
      if(norm(aij=A(i,j))>tol2)
       {
         cl[nbcoeff]=j;
         a[nbcoeff]=aij;
         nbcoeff++;
       }
    }
   lg[this->n]=nbcoeff;

  
}
template <class R> 
 MatriceMorse<R>::MatriceMorse(const int  nn,const R *aa)
    :MatriceCreuse<R>(nn),solver(0) 
      {
  symetrique = true;
  this->dummy=false;
  this->n=nn;
  nbcoef=this->n;
  a=new R[this->n] ;
  lg=new int [this->n+1];
  cl=new int [this->n];
  for(int i=0;i<this->n;i++)
   {
    lg[i]=i;
    cl[i]=i;
    a[i]=aa[i];      
      }
lg[this->n]=this->n;
}

template<class R>
template<class K>
 MatriceMorse<R>::MatriceMorse(const MatriceMorse<K> & A)
   : MatriceCreuse<R>(A.n,A.m,A.dummy),nbcoef(A.nbcoef),      
     symetrique(A.symetrique),       
     a(new R[nbcoef]),
     lg(new int [this->n+1]),
     cl(new int[nbcoef]),
     solver(0)
{
  ffassert(a && lg &&  cl);
  for (int i=0;i<=this->n;i++)
    lg[i]=A.lg[i];
  for (int k=0;k<nbcoef;k++)
    {
      cl[k]=A.cl[k];
      a[k]=A.a[k];
    }
  
}



template <class R> 
int MatriceMorse<R>::size() const 
{
  return nbcoef*(sizeof(int)+sizeof(R))+ sizeof(int)*(this->n+1);
}

inline int WhichMatrix(istream & f)
{
    string line;
    while ( isspace(f.peek()))
	f.get(); 
    if  ( f.peek() =='#' )
    {
	line="";
	while ( f.good()  )
	{
	    char c=f.get();
	    if(c=='\n' || c=='\r') { break;}
	    line += c;
	}
	if( line.find("(Morse)")) 
	    return 2; // morse 
	else 
	    return 0; 
    }
  return 0;   
}
template <class R>
  MatriceMorse<R>::MatriceMorse(istream & f)
:  MatriceCreuse<R>(0,0,0),nbcoef(0),
a(0),
lg(0),
cl(0),

solver(0)
{
      string line;
      int k=0;
      while ( isspace(f.peek()))
	      f.get(); 
      while ( f.peek() =='#' )
      {
	line="";
	while ( f.good()  )
	{
	    char c=f.get();
	    if(c=='\n' || c=='\r') { break;}
	    line += c;
	}
	if( f.peek()=='\n' || f.peek()=='\r') f.get();
	if(verbosity>9)
	 cout << "Read matrice: "<< k << " :"   << line << endl;
	k++;    
      }
      
      f >> this->n >> this->m >> symetrique >>nbcoef;
      if(verbosity>3)
      cout << " read mat: " <<  this->n << " " <<  this->m << " " << symetrique << " " << nbcoef <<endl;
      lg= new int [this->n+1];
      cl= new int[nbcoef];
      a= new R[nbcoef];
      ffassert(f.good() && lg && a && cl );
      int i,j,i0,j0;
      i0=-1;j0=2000000000;
      R aij;
      int imx=-2000000000, jmx=-2000000000;
      int imn= 2000000000, jmn= 2000000000;
      
      for (int k =0;k<nbcoef; ++k)
      {
	  f >> i >> j >> aij;
	  ffassert(f.good() );
	  i--;j--;
	  imx=max(imx,i);
	  imx=max(jmx,j);
	  imn=min(imn,i);
	  imn=min(jmn,j);
	  //cout << i << " " << j << " " << aij << endl;
	  if(i0!=i) {j0=-1;lg[i]=k;}
	  ffassert(i0<=i && j0<j);
	  lg[i+1]=k+1;
	  cl[k]=j;
	  a[k]=aij;
	  j0=j;i0=i;
      }
      ffassert( imx < this->n && jmx < this->m );
      ffassert( imn >=0 && jmn >=0);
      
}

template <class R> 
ostream& MatriceMorse<R>::dump(ostream & f) const 
{