op_dsdphi.C
/*
* Copyright (c) 1999-2000 Jean-Alain Marck
* Copyright (c) 1999-2001 Eric Gourgoulhon
*
* This file is part of LORENE.
*
* LORENE 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.
*
* LORENE 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 LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char op_dsdphi_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_dsdphi.C,v 1.6 2014/10/13 08:53:25 j_novak Exp $" ;
/*
* Ensemble des routines de base pour la derivation par rapport a phi
* (Utilisation interne)
*
* void _dsdphi_XXXX(Tbl * t, int & b)
* t pointeur sur le Tbl d'entree-sortie
* b base spectrale
*
*/
/*
* $Id: op_dsdphi.C,v 1.6 2014/10/13 08:53:25 j_novak Exp $
* $Log: op_dsdphi.C,v $
* Revision 1.6 2014/10/13 08:53:25 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.5 2013/04/25 15:46:06 j_novak
* Added special treatment in the case np = 1, for type_p = NONSYM.
*
* Revision 1.4 2006/03/10 12:45:38 j_novak
* Use of C++-style cast.
*
* Revision 1.3 2003/12/19 16:21:48 j_novak
* Shadow hunt
*
* Revision 1.2 2002/09/09 13:00:40 e_gourgoulhon
* Modification of declaration of Fortran 77 prototypes for
* a better portability (in particular on IBM AIX systems):
* All Fortran subroutine names are now written F77_* and are
* defined in the new file C++/Include/proto_f77.h.
*
* Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
* LORENE
*
* Revision 2.8 2000/11/14 15:08:16 eric
* Traitement du cas np=1 dans P_COSSIN_I.
*
* Revision 2.7 2000/10/04 15:54:09 eric
* *** empty log message ***
*
* Revision 2.6 2000/10/04 12:50:38 eric
* Ajout de la base P_COSSIN_I.
*
* Revision 2.5 2000/10/04 11:50:32 eric
* Ajout d' abort() dans le cas non prevu.
*
* Revision 2.4 2000/08/18 13:20:01 eric
* Traitement du cas np = 1.
*
* Revision 2.3 2000/02/24 16:41:03 eric
* Initialisation a zero du tableau xo avant le calcul.
*
* Revision 2.2 1999/11/15 16:38:03 eric
* Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
*
* Revision 2.1 1999/02/23 11:36:34 hyc
* *** empty log message ***
*
* Revision 2.0 1999/02/22 17:24:59 hyc
* *** empty log message ***
*
* $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_dsdphi.C,v 1.6 2014/10/13 08:53:25 j_novak Exp $
*
*/
// Headers Lorene
#include "tbl.h"
#include "proto_f77.h"
// Routine pour les cas non prevus
//--------------------------------
namespace Lorene {
void _dsdphi_pas_prevu(Tbl* , int & b) {
cout << "Unknown phi basis in Mtbl_cf::dsdp() !" << endl ;
cout << " basis: " << hex << b << endl ;
abort() ;
}
//------------------//
// cas P_COSSIN //
//------------------//
void _dsdphi_p_cossin(Tbl* tb, int & )
{
// Peut-etre rien a faire ?
if (tb->get_etat() == ETATZERO) {
return ;
}
// Protection
assert(tb->get_etat() == ETATQCQ) ;
// Pour le confort
int nr = (tb->dim).dim[0] ; // Nombre
int nt = (tb->dim).dim[1] ; // de points
int np = (tb->dim).dim[2] ; // physiques REELS
np = np - 2 ; // Nombre de points physiques
// Cas particulier de la symetrie axiale :
if (np == 1) {
tb->set_etat_zero() ;
return ;
}
// Variables statiques
static double* cx = 0 ;
static int np_pre =0 ;
// Test sur np pour initialisation eventuelle
//#pragma critical (loch_dsdphi_p_cossin)
{
if (np > np_pre) {
np_pre = np ;
cx = reinterpret_cast<double*>(realloc(cx, (np+2) * sizeof(double))) ;
for (int i=0 ; i<np+2 ; i += 2) {
cx[i] = - (i/2) ;
cx[i+1] = (i/2) ;
}
}
} // Fin de region critique
// pt. sur le tableau de double resultat
double* xo = new double[(tb->dim).taille] ;
// Initialisation a zero :
for (int i=0; i<(tb->dim).taille; i++) {
xo[i] = 0 ;
}
// On y va...
double* xi = tb->t ;
double* xci = xi ; // Pointeurs
double* xco = xo ; // courants
int k, j ;
// k = 0: inutile, deviendra k=1
xci += nr*nt ; // Positionnement
xco += nr*nt ;
// k = 1: 0
for (j=0 ; j<nr*nt ; j++) {
xco[j] = 0 ;
} // Fin de la boucle sur r * theta
xci += nr*nt ; // Positionnement
xco += nr*nt ;
// 2 =< k <= np-1
for (k=2 ; k<np ; k++) {
for (j=0 ; j<nr*nt ; j++) {
xco[j] = cx[k] * xci[j] ;
} // Fin de la boucle sur r * theta
xci += nr*nt ; // Positionnement
xco += nr*nt ;
} // Fin de la boucle sur phi
// k = np: inutile, deviendra k=np+1
xci += nr*nt ; // Positionnement
xco += nr*nt ;
// k = np+1
for (j=0 ; j<nr*nt ; j++) {
xco[j] = 0 ;
} // Fin de la boucle sur r * theta
// Inversion des sinus et cosinus (appel a dswap de BLAS)
int nbr = nr*nt ;
int inc = 1 ;
double* p1 = xo ;
double* p2 = p1 + nr*nt ;
for (k=0 ; k<np+1 ; k +=2, p1 += 2*nr*nt, p2 += 2*nr*nt) {
F77_dswap(&nbr, p1, &inc, p2, &inc) ;
}
// On remet les choses la ou il faut
delete [] tb->t ;
tb->t = xo ;
// base de developpement
// inchangee
}
//------------------//
// cas P_COSSIN_P //
//------------------//
void _dsdphi_p_cossin_p(Tbl* tb, int & )
{
// Peut-etre rien a faire ?
if (tb->get_etat() == ETATZERO) {
return ;
}
// Protection
assert(tb->get_etat() == ETATQCQ) ;
// Pour le confort
int nr = (tb->dim).dim[0] ; // Nombre
int nt = (tb->dim).dim[1] ; // de points
int np = (tb->dim).dim[2] ; // physiques REELS
np = np - 2 ; // Nombre de points physiques
// Cas particulier de la symetrie axiale :
if (np == 1) {
tb->set_etat_zero() ;
return ;
}
// Variables statiques
static double* cx = 0 ;
static int np_pre =0 ;
// Test sur np pour initialisation eventuelle
//#pragma critical (loch_dsdphi_p_cossin_p)
{
if (np > np_pre) {
np_pre = np ;
cx = reinterpret_cast<double*>(realloc(cx, (np+2) * sizeof(double))) ;
int i ;
for (i=0 ; i<np+2 ; i += 2) {
cx[i] = - (i/2) ;
cx[i+1] = (i/2) ;
}
for (i=0 ; i<np+2 ; i++) {
cx[i] = 2 * cx[i] ;
}
}
} // Fin de region critique
// pt. sur le tableau de double resultat
double* xo = new double[(tb->dim).taille] ;
// Initialisation a zero :
for (int i=0; i<(tb->dim).taille; i++) {
xo[i] = 0 ;
}
// On y va...
double* xi = tb->t ;
double* xci = xi ; // Pointeurs
double* xco = xo ; // courants
int k, j ;
// k = 0: inutile, deviendra k=1
xci += nr*nt ; // Positionnement
xco += nr*nt ;
// k = 1: 0
for (j=0 ; j<nr*nt ; j++) {
xco[j] = 0 ;
} // Fin de la boucle sur r * theta
xci += nr*nt ; // Positionnement
xco += nr*nt ;
// 2 =< k <= np-1
for (k=2 ; k<np ; k++) {
for (j=0 ; j<nr*nt ; j++) {
xco[j] = cx[k] * xci[j] ;
} // Fin de la boucle sur r * theta
xci += nr*nt ; // Positionnement
xco += nr*nt ;
} // Fin de la boucle sur phi
// k = np: inutile, deviendra k=np+1
xci += nr*nt ; // Positionnement
xco += nr*nt ;
// k = np+1
for (j=0 ; j<nr*nt ; j++) {
xco[j] = 0 ;
} // Fin de la boucle sur r * theta
// Inversion des sinus et cosinus (appel a dswap de BLAS)
int nbr = nr*nt ;
int inc = 1 ;
double* p1 = xo ;
double* p2 = p1 + nr*nt ;
for (k=0 ; k<np+1 ; k +=2, p1 += 2*nr*nt, p2 += 2*nr*nt) {
F77_dswap(&nbr, p1, &inc, p2, &inc) ;
}
// On remet les choses la ou il faut
delete [] tb->t ;
tb->t = xo ;
// base de developpement
// inchangee
}
//------------------//
// cas P_COSSIN_I //
//------------------//
void _dsdphi_p_cossin_i(Tbl* tb, int & )
{
// Peut-etre rien a faire ?
if (tb->get_etat() == ETATZERO) {
return ;
}
// Protection
assert(tb->get_etat() == ETATQCQ) ;
// Pour le confort
int nr = (tb->dim).dim[0] ; // Nombre
int nt = (tb->dim).dim[1] ; // de points
int np = (tb->dim).dim[2] ; // physiques REELS
np = np - 2 ; // Nombre de points physiques
int ntnr = nt*nr ; // increment en phi
// Cas particulier de la symetrie axiale :
if (np == 1) {
tb->set_etat_zero() ;
return ;
}
// pt. sur le tableau de double resultat
double* xo = new double[(tb->dim).taille] ;
// pt. sur le tableau de double entree
const double* xi = tb->t ;
// k = 0 : resultat sur cos(phi)
// ------------------------------
double* xco = xo ; // coef de cos(phi)
const double* xci = xi + 2*ntnr ; // coef de sin(phi)
for (int i=0; i<ntnr; i++) {
xco[i] = xci[i] ;
}
// k = 1 : mise a zero
// --------------------
xco = xo + ntnr ;
for (int i=0; i<ntnr; i++) {
xco[i] = 0 ;
}
// k = 2 : resultat sur sin(phi)
// ------------------------------
xco = xo + 2*ntnr ; // coef de sin(phi)
xci = xi ; // coef de cos(phi)
for (int i=0; i<ntnr; i++) {
xco[i] = - xci[i] ;
}
// k >= 3
// ------
for (int k=3; k<np; k+=2) {
// resultat sur cos(k phi)
// -----------------------
xco = xo + k*ntnr ; // coef de cos(k phi)
xci = xi + (k+1)*ntnr ; // coef de sin(k phi)
for (int i=0; i<ntnr; i++) {
xco[i] = k * xci[i] ;
}
// resultat sur sin(k phi)
// -----------------------
xco = xo + (k+1)*ntnr ; // coef de sin(k phi)
xci = xi + k*ntnr ; // coef de cos(k phi)
for (int i=0; i<ntnr; i++) {
xco[i] = - k * xci[i] ;
}
}
// k = np+1 : mise a zero
// -----------------------
xco = xo + (np+1)*ntnr ;
for (int i=0; i<ntnr; i++) {
xco[i] = 0 ;
}
// On remet les choses la ou il faut
// ---------------------------------
delete [] tb->t ;
tb->t = xo ;
// base de developpement
// inchangee
}
}
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene
Lorene prototypes.
Definition app_hor.h:64