OpenCascade/occt
1
0
Fork 0
mirror of https://git.dev.opencascade.org/repos/occt.git synced 2025-04-21 10:13:43 +03:00
Code
cfdd99315c
occt/src/AdvApp2Var/AdvApp2Var_MathBase.cxx
abv ed4415982c 0024624: Lost word in license statement in source files 
License statement text corrected; compiler warnings caused by Bison 2.41 disabled for MSVC; a few other compiler warnings on 54-bit Windows eliminated by appropriate type cast
Wrong license statements corrected in several files.
Copyright and license statements added in XSD and GLSL files.
Copyright year updated in some files.
Obsolete documentation files removed from DrawResources.
2014-03-25 16:38:55 +04:00

11301 lines
292 KiB
C++
Raw Blame History

// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.
// AdvApp2Var_MathBase.cxx
#include <math.h>
#include <AdvApp2Var_SysBase.hxx>
#include <AdvApp2Var_Data_f2c.hxx>
#include <AdvApp2Var_MathBase.hxx>
#include <AdvApp2Var_Data.hxx>
// statics
static
int mmchole_(integer *mxcoef,
integer *dimens,
doublereal *amatri,
integer *aposit,
integer *posuiv,
doublereal *chomat,
integer *iercod);
static
int mmrslss_(integer *mxcoef,
integer *dimens,
doublereal *smatri,
integer *sposit,
integer *posuiv,
doublereal *mscnmbr,
doublereal *soluti,
integer *iercod);
static
int mfac_(doublereal *f,
integer *n);
static
int mmaper0_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *crvlgd,
integer *ncfnew,
doublereal *ycvmax,
doublereal *errmax);
static
int mmaper2_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *crvjac,
integer *ncfnew,
doublereal *ycvmax,
doublereal *errmax);
static
int mmaper4_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *crvjac,
integer *ncfnew,
doublereal *ycvmax,
doublereal *errmax);
static
int mmaper6_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *crvjac,
integer *ncfnew,
doublereal *ycvmax,
doublereal *errmax);
static
int mmarc41_(integer *ndimax,
integer *ndimen,
integer *ncoeff,
doublereal *crvold,
doublereal *upara0,
doublereal *upara1,
doublereal *crvnew,
integer *iercod);
static
int mmatvec_(integer *nligne,
integer *ncolon,
integer *gposit,
integer *gnstoc,
doublereal *gmatri,
doublereal *vecin,
integer *deblig,
doublereal *vecout,
integer *iercod);
static
int mmcvstd_(integer *ncofmx,
integer *ndimax,
integer *ncoeff,
integer *ndimen,
doublereal *crvcan,
doublereal *courbe);
static
int mmdrvcb_(integer *ideriv,
integer *ndim,
integer *ncoeff,
doublereal *courbe,
doublereal *tparam,
doublereal *tabpnt,
integer *iercod);
static
int mmexthi_(integer *ndegre,
doublereal *hwgaus);
static
int mmextrl_(integer *ndegre,
doublereal *rootlg);
static
int mmherm0_(doublereal *debfin,
integer *iercod);
static
int mmherm1_(doublereal *debfin,
integer *ordrmx,
integer *iordre,
doublereal *hermit,
integer *iercod);
static
int mmloncv_(integer *ndimax,
integer *ndimen,
integer *ncoeff,
doublereal *courbe,
doublereal *tdebut,
doublereal *tfinal,
doublereal *xlongc,
integer *iercod);
static
int mmpojac_(doublereal *tparam,
integer *iordre,
integer *ncoeff,
integer *nderiv,
doublereal *valjac,
integer *iercod);
static
int mmrslw_(integer *normax,
integer *nordre,
integer *ndimen,
doublereal *epspiv,
doublereal *abmatr,
doublereal *xmatri,
integer *iercod);
static
int mmtmave_(integer *nligne,
integer *ncolon,
integer *gposit,
integer *gnstoc,
doublereal *gmatri,
doublereal *vecin,
doublereal *vecout,
integer *iercod);
static
int mmtrpj0_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *epsi3d,
doublereal *crvlgd,
doublereal *ycvmax,
doublereal *epstrc,
integer *ncfnew);
static
int mmtrpj2_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *epsi3d,
doublereal *crvlgd,
doublereal *ycvmax,
doublereal *epstrc,
integer *ncfnew);
static
int mmtrpj4_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *epsi3d,
doublereal *crvlgd,
doublereal *ycvmax,
doublereal *epstrc,
integer *ncfnew);
static
int mmtrpj6_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *epsi3d,
doublereal *crvlgd,
doublereal *ycvmax,
doublereal *epstrc,
integer *ncfnew);
static
integer pow__ii(integer *x,
integer *n);
static
int mvcvin2_(integer *ncoeff,
doublereal *crvold,
doublereal *crvnew,
integer *iercod);
static
int mvcvinv_(integer *ncoeff,
doublereal *crvold,
doublereal *crvnew,
integer *iercod);
static
int mvgaus0_(integer *kindic,
doublereal *urootl,
doublereal *hiltab,
integer *nbrval,
integer *iercod);
static
int mvpscr2_(integer *ncoeff,
doublereal *curve2,
doublereal *tparam,
doublereal *pntcrb);
static
int mvpscr3_(integer *ncoeff,
doublereal *curve2,
doublereal *tparam,
doublereal *pntcrb);
static struct {
doublereal eps1, eps2, eps3, eps4;
integer niterm, niterr;
} mmprcsn_;
static struct {
doublereal tdebut, tfinal, verifi, cmherm[576];
} mmcmher_;
//=======================================================================
//function : AdvApp2Var_MathBase::mdsptpt_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mdsptpt_(integer *ndimen,
doublereal *point1,
doublereal *point2,
doublereal *distan)
{
integer c__8 = 8;
/* System generated locals */
integer i__1;
doublereal d__1;
/* Local variables */
integer i__;
doublereal* differ = 0;
integer ier;
intptr_t iofset, j;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* CALCULATE DISTANCE BETWEEN TWO POINTS */
/* KEYWORDS : */
/* ----------- */
/* DISTANCE,POINT. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMEN: Space Dimension. */
/* POINT1: Table of coordinates of the 1st point. */
/* POINT2: Table of coordinates of the 2nd point. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* DISTAN: Distance between 2 points. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* **********************************************************************
*/
/* ***********************************************************************
*/
/* INITIALISATION */
/* ***********************************************************************
*/
/* Parameter adjustment */
--point2;
--point1;
/* Function Body */
iofset = 0;
ier = 0;
/* ***********************************************************************
*/
/* TRAITEMENT */
/* ***********************************************************************
*/
AdvApp2Var_SysBase anAdvApp2Var_SysBase;
if (*ndimen > 100) {
anAdvApp2Var_SysBase.mcrrqst_(&c__8, ndimen, differ, &iofset, &ier);
}
/* --- If allocation is refused, the trivial method is applied. */
if (ier > 0) {
*distan = 0.;
i__1 = *ndimen;
for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing 2nd power */
d__1 = point1[i__] - point2[i__];
*distan += d__1 * d__1;
}
*distan = sqrt(*distan);
/* --- Otherwise MZSNORM is used to minimize the risks of overflow
*/
} else {
i__1 = *ndimen;
for (i__ = 1; i__ <= i__1; ++i__) {
j=iofset + i__ - 1;
differ[j] = point2[i__] - point1[i__];
}
*distan = AdvApp2Var_MathBase::mzsnorm_(ndimen, &differ[iofset]);
}
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
/* --- Dynamic Desallocation */
if (iofset != 0) {
anAdvApp2Var_SysBase.mcrdelt_(&c__8, ndimen, differ, &iofset, &ier);
}
return 0 ;
} /* mdsptpt_ */
//=======================================================================
//function : mfac_
//purpose :
//=======================================================================
int mfac_(doublereal *f,
integer *n)
{
/* System generated locals */
integer i__1;
/* Local variables */
integer i__;
/* FORTRAN CONFORME AU TEXT */
/* CALCUL DE MFACTORIEL N */
/* Parameter adjustments */
--f;
/* Function Body */
f[1] = (float)1.;
i__1 = *n;
for (i__ = 2; i__ <= i__1; ++i__) {
/* L10: */
f[i__] = i__ * f[i__ - 1];
}
return 0;
} /* mfac_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmapcmp_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmapcmp_(integer *ndim,
integer *ncofmx,
integer *ncoeff,
doublereal *crvold,
doublereal *crvnew)
{
/* System generated locals */
integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
i__2;
/* Local variables */
integer ipair, nd, ndegre, impair, ibb, idg;
//extern int mgsomsg_();//mgenmsg_(),
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Compression of curve CRVOLD in a table of */
/* coeff. of even : CRVNEW(*,0,*) */
/* and uneven range : CRVNEW(*,1,*). */
/* KEYWORDS : */
/* ----------- */
/* COMPRESSION,CURVE. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIM : Space Dimension. */
/* NCOFMX : Max nb of coeff. of the curve to compress. */
/* NCOEFF : Max nb of coeff. of the compressed curve. */
/* CRVOLD : The curve (0:NCOFMX-1,NDIM) to compress. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* CRVNEW : Curve compacted in (0:(NCOEFF-1)/2,0,NDIM) (containing
*/
/* even terms) and in (0:(NCOEFF-1)/2,1,NDIM) */
/* (containing uneven terms). */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* This routine is useful to prepare coefficients of a */
/* curve in an orthogonal base (Legendre or Jacobi) before */
/* calculating the coefficients in the canonical; base [-1,1] by */
/* MMJACAN. */
/* ***********************************************************************
*/
/* Name of the routine */
/* Parameter adjustments */
crvold_dim1 = *ncofmx;
crvold_offset = crvold_dim1;
crvold -= crvold_offset;
crvnew_dim1 = (*ncoeff - 1) / 2 + 1;
crvnew_offset = crvnew_dim1 << 1;
crvnew -= crvnew_offset;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 3) {
AdvApp2Var_SysBase::mgenmsg_("MMAPCMP", 7L);
}
ndegre = *ncoeff - 1;
i__1 = *ndim;
for (nd = 1; nd <= i__1; ++nd) {
ipair = 0;
i__2 = ndegre / 2;
for (idg = 0; idg <= i__2; ++idg) {
crvnew[idg + (nd << 1) * crvnew_dim1] = crvold[ipair + nd *
crvold_dim1];
ipair += 2;
/* L200: */
}
if (ndegre < 1) {
goto L400;
}
impair = 1;
i__2 = (ndegre - 1) / 2;
for (idg = 0; idg <= i__2; ++idg) {
crvnew[idg + ((nd << 1) + 1) * crvnew_dim1] = crvold[impair + nd *
crvold_dim1];
impair += 2;
/* L300: */
}
L400:
/* L100: */
;
}
/* ---------------------------------- The end ---------------------------
*/
if (ibb >= 3) {
AdvApp2Var_SysBase::mgsomsg_("MMAPCMP", 7L);
}
return 0;
} /* mmapcmp_ */
//=======================================================================
//function : mmaper0_
//purpose :
//=======================================================================
int mmaper0_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *crvlgd,
integer *ncfnew,
doublereal *ycvmax,
doublereal *errmax)
{
/* System generated locals */
integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer ncut;
doublereal bidon;
integer ii, nd;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Calculate the max error of approximation done when */
/* only the first NCFNEW coefficients of a curve are preserved.
*/
/* Degree NCOEFF-1 written in the base of Legendre (Jacobi */
/* of order 0). */
/* KEYWORDS : */
/* ----------- */
/* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX : Max. degree of the curve. */
/* NDIMEN : Space dimension. */
/* NCOEFF : Degree +1 of the curve. */
/* CRVLGD : Curve the degree which of should be lowered. */
/* NCFNEW : Degree +1 of the resulting polynom. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* YCVMAX : Auxiliary Table (max error on each dimension).
*/
/* ERRMAX : Precision of the approximation. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ***********************************************************************
*/
/* ------------------- Init to calculate an error -----------------------
*/
/* Parameter adjustments */
--ycvmax;
crvlgd_dim1 = *ncofmx;
crvlgd_offset = crvlgd_dim1 + 1;
crvlgd -= crvlgd_offset;
/* Function Body */
i__1 = *ndimen;
for (ii = 1; ii <= i__1; ++ii) {
ycvmax[ii] = 0.;
/* L100: */
}
/* ------ Minimum that can be reached : Stop at 1 or NCFNEW ------
*/
ncut = 1;
if (*ncfnew + 1 > ncut) {
ncut = *ncfnew + 1;
}
/* -------------- Elimination of high degree coefficients-----------
*/
/* ----------- Loop on the series of Legendre: NCUT --> NCOEFF --------
*/
i__1 = *ncoeff;
for (ii = ncut; ii <= i__1; ++ii) {
/* Factor of renormalization (Maximum of Li(t)). */
bidon = ((ii - 1) * 2. + 1.) / 2.;
bidon = sqrt(bidon);
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
ycvmax[nd] += (d__1 = crvlgd[ii + nd * crvlgd_dim1], advapp_abs(d__1)) *
bidon;
/* L310: */
}
/* L300: */
}
/* -------------- The error is the norm of the vector error ---------------
*/
*errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
/* --------------------------------- Fin --------------------------------
*/
return 0;
} /* mmaper0_ */
//=======================================================================
//function : mmaper2_
//purpose :
//=======================================================================
int mmaper2_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *crvjac,
integer *ncfnew,
doublereal *ycvmax,
doublereal *errmax)
{
/* Initialized data */
static doublereal xmaxj[57] = { .9682458365518542212948163499456,
.986013297183269340427888048593603,
1.07810420343739860362585159028115,
1.17325804490920057010925920756025,
1.26476561266905634732910520370741,
1.35169950227289626684434056681946,
1.43424378958284137759129885012494,
1.51281316274895465689402798226634,
1.5878364329591908800533936587012,
1.65970112228228167018443636171226,
1.72874345388622461848433443013543,
1.7952515611463877544077632304216,
1.85947199025328260370244491818047,
1.92161634324190018916351663207101,
1.98186713586472025397859895825157,
2.04038269834980146276967984252188,
2.09730119173852573441223706382076,
2.15274387655763462685970799663412,
2.20681777186342079455059961912859,
2.25961782459354604684402726624239,
2.31122868752403808176824020121524,
2.36172618435386566570998793688131,
2.41117852396114589446497298177554,
2.45964731268663657873849811095449,
2.50718840313973523778244737914028,
2.55385260994795361951813645784034,
2.59968631659221867834697883938297,
2.64473199258285846332860663371298,
2.68902863641518586789566216064557,
2.73261215675199397407027673053895,
2.77551570192374483822124304745691,
2.8177699459714315371037628127545,
2.85940333797200948896046563785957,
2.90044232019793636101516293333324,
2.94091151970640874812265419871976,
2.98083391718088702956696303389061,
3.02023099621926980436221568258656,
3.05912287574998661724731962377847,
3.09752842783622025614245706196447,
3.13546538278134559341444834866301,
3.17295042316122606504398054547289,
3.2099992681699613513775259670214,
3.24662674946606137764916854570219,
3.28284687953866689817670991319787,
3.31867291347259485044591136879087,
3.35411740487202127264475726990106,
3.38919225660177218727305224515862,
3.42390876691942143189170489271753,
3.45827767149820230182596660024454,
3.49230918177808483937957161007792,
3.5260130200285724149540352829756,
3.55939845146044235497103883695448,
3.59247431368364585025958062194665,
3.62524904377393592090180712976368,
3.65773070318071087226169680450936,
3.68992700068237648299565823810245,
3.72184531357268220291630708234186 };
/* System generated locals */
integer crvjac_dim1, crvjac_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer idec, ncut;
doublereal bidon;
integer ii, nd;
/* ***********************************************************************
*/
/* FONCTION : */
/* ---------- */
/* Calculate max approximation error i faite lorsque l' on */
/* ne conserve que les premiers NCFNEW coefficients d' une courbe
*/
/* de degre NCOEFF-1 ecrite dans la base de Jacobi d' ordre 2. */
/* KEYWORDS : */
/* ----------- */
/* JACOBI, POLYGON, APPROXIMATION, ERROR. */
/**/
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX : Max. degree of the curve. */
/* NDIMEN : Space dimension. */
/* NCOEFF : Degree +1 of the curve. */
/* CRVLGD : Curve the degree which of should be lowered. */
/* NCFNEW : Degree +1 of the resulting polynom. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* YCVMAX : Auxiliary Table (max error on each dimension).
*/
/* ERRMAX : Precision of the approximation. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ------------------ Table of maximums of (1-t2)*Ji(t) ----------------
*/
/* Parameter adjustments */
--ycvmax;
crvjac_dim1 = *ncofmx;
crvjac_offset = crvjac_dim1 + 1;
crvjac -= crvjac_offset;
/* Function Body */
/* ------------------- Init for error calculation -----------------------
*/
i__1 = *ndimen;
for (ii = 1; ii <= i__1; ++ii) {
ycvmax[ii] = 0.;
/* L100: */
}
/* ------ Min. Degree that can be attained : Stop at 3 or NCFNEW ------
*/
idec = 3;
/* Computing MAX */
i__1 = idec, i__2 = *ncfnew + 1;
ncut = advapp_max(i__1,i__2);
/* -------------- Removal of coefficients of high degree -----------
*/
/* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
*/
i__1 = *ncoeff;
for (ii = ncut; ii <= i__1; ++ii) {
/* Factor of renormalization. */
bidon = xmaxj[ii - idec];
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
bidon;
/* L310: */
}
/* L300: */
}
/* -------------- The error is the norm of the vector error ---------------
*/
*errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
/* --------------------------------- Fin --------------------------------
*/
return 0;
} /* mmaper2_ */
/* MAPER4.f -- translated by f2c (version 19960827).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
/* Subroutine */
//=======================================================================
//function : mmaper4_
//purpose :
//=======================================================================
int mmaper4_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *crvjac,
integer *ncfnew,
doublereal *ycvmax,
doublereal *errmax)
{
/* Initialized data */
static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
1.05299572648705464724876659688996,
1.0949715351434178709281698645813,
1.15078388379719068145021100764647,
1.2094863084718701596278219811869,
1.26806623151369531323304177532868,
1.32549784426476978866302826176202,
1.38142537365039019558329304432581,
1.43575531950773585146867625840552,
1.48850442653629641402403231015299,
1.53973611681876234549146350844736,
1.58953193485272191557448229046492,
1.63797820416306624705258190017418,
1.68515974143594899185621942934906,
1.73115699602477936547107755854868,
1.77604489805513552087086912113251,
1.81989256661534438347398400420601,
1.86276344480103110090865609776681,
1.90471563564740808542244678597105,
1.94580231994751044968731427898046,
1.98607219357764450634552790950067,
2.02556989246317857340333585562678,
2.06433638992049685189059517340452,
2.10240936014742726236706004607473,
2.13982350649113222745523925190532,
2.17661085564771614285379929798896,
2.21280102016879766322589373557048,
2.2484214321456956597803794333791,
2.28349755104077956674135810027654,
2.31805304852593774867640120860446,
2.35210997297725685169643559615022,
2.38568889602346315560143377261814,
2.41880904328694215730192284109322,
2.45148841120796359750021227795539,
2.48374387161372199992570528025315,
2.5155912654873773953959098501893,
2.54704548720896557684101746505398,
2.57812056037881628390134077704127,
2.60882970619319538196517982945269,
2.63918540521920497868347679257107,
2.66919945330942891495458446613851,
2.69888301230439621709803756505788,
2.72824665609081486737132853370048,
2.75730041251405791603760003778285,
2.78605380158311346185098508516203,
2.81451587035387403267676338931454,
2.84269522483114290814009184272637,
2.87060005919012917988363332454033,
2.89823818258367657739520912946934,
2.92561704377132528239806135133273,
2.95274375377994262301217318010209,
2.97962510678256471794289060402033,
3.00626759936182712291041810228171,
3.03267744830655121818899164295959,
3.05886060707437081434964933864149 };
/* System generated locals */
integer crvjac_dim1, crvjac_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer idec, ncut;
doublereal bidon;
integer ii, nd;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Calculate the max. error of approximation made when */
/* only first NCFNEW coefficients of a curve are preserved
*/
/* degree NCOEFF-1 is written in the base of Jacobi of order 4. */
/* KEYWORDS : */
/* ----------- */
/* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX : Max. degree of the curve. */
/* NDIMEN : Space dimension. */
/* NCOEFF : Degree +1 of the curve. */
/* CRVJAC : Curve the degree which of should be lowered. */
/* NCFNEW : Degree +1 of the resulting polynom. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* YCVMAX : Auxiliary Table (max error on each dimension).
*/
/* ERRMAX : Precision of the approximation. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ***********************************************************************
*/
/* ---------------- Table of maximums of ((1-t2)2)*Ji(t) ---------------
*/
/* Parameter adjustments */
--ycvmax;
crvjac_dim1 = *ncofmx;
crvjac_offset = crvjac_dim1 + 1;
crvjac -= crvjac_offset;
/* Function Body */
/* ------------------- Init for error calculation -----------------------
*/
i__1 = *ndimen;
for (ii = 1; ii <= i__1; ++ii) {
ycvmax[ii] = 0.;
/* L100: */
}
/* ------ Min. Degree that can be attained : Stop at 5 or NCFNEW ------
*/
idec = 5;
/* Computing MAX */
i__1 = idec, i__2 = *ncfnew + 1;
ncut = advapp_max(i__1,i__2);
/* -------------- Removal of high degree coefficients -----------
*/
/* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
*/
i__1 = *ncoeff;
for (ii = ncut; ii <= i__1; ++ii) {
/* Factor of renormalisation. */
bidon = xmaxj[ii - idec];
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
bidon;
/* L310: */
}
/* L300: */
}
/* -------------- The error is the norm of the error vector ---------------
*/
*errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
/* --------------------------------- End --------------------------------
*/
return 0;
} /* mmaper4_ */
//=======================================================================
//function : mmaper6_
//purpose :
//=======================================================================
int mmaper6_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *crvjac,
integer *ncfnew,
doublereal *ycvmax,
doublereal *errmax)
{
/* Initialized data */
static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
1.11626917091567929907256116528817,
1.1327140810290884106278510474203,
1.1679452722668028753522098022171,
1.20910611986279066645602153641334,
1.25228283758701572089625983127043,
1.29591971597287895911380446311508,
1.3393138157481884258308028584917,
1.3821288728999671920677617491385,
1.42420414683357356104823573391816,
1.46546895108549501306970087318319,
1.50590085198398789708599726315869,
1.54550385142820987194251585145013,
1.58429644271680300005206185490937,
1.62230484071440103826322971668038,
1.65955905239130512405565733793667,
1.69609056468292429853775667485212,
1.73193098017228915881592458573809,
1.7671112206990325429863426635397,
1.80166107681586964987277458875667,
1.83560897003644959204940535551721,
1.86898184653271388435058371983316,
1.90180515174518670797686768515502,
1.93410285411785808749237200054739,
1.96589749778987993293150856865539,
1.99721027139062501070081653790635,
2.02806108474738744005306947877164,
2.05846864831762572089033752595401,
2.08845055210580131460156962214748,
2.11802334209486194329576724042253,
2.14720259305166593214642386780469,
2.17600297710595096918495785742803,
2.20443832785205516555772788192013,
2.2325216999457379530416998244706,
2.2602654243075083168599953074345,
2.28768115912702794202525264301585,
2.3147799369092684021274946755348,
2.34157220782483457076721300512406,
2.36806787963276257263034969490066,
2.39427635443992520016789041085844,
2.42020656255081863955040620243062,
2.44586699364757383088888037359254,
2.47126572552427660024678584642791,
2.49641045058324178349347438430311,
2.52130850028451113942299097584818,
2.54596686772399937214920135190177,
2.5703922285006754089328998222275,
2.59459096001908861492582631591134,
2.61856915936049852435394597597773,
2.64233265984385295286445444361827,
2.66588704638685848486056711408168,
2.68923766976735295746679957665724,
2.71238965987606292679677228666411 };
/* System generated locals */
integer crvjac_dim1, crvjac_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer idec, ncut;
doublereal bidon;
integer ii, nd;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Calculate the max. error of approximation made when */
/* only first NCFNEW coefficients of a curve are preserved
*/
/* degree NCOEFF-1 is written in the base of Jacobi of order 6. */
/* KEYWORDS : */
/* ----------- */
/* JACOBI,POLYGON,APPROXIMATION,ERROR. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX : Max. degree of the curve. */
/* NDIMEN : Space dimension. */
/* NCOEFF : Degree +1 of the curve. */
/* CRVJAC : Curve the degree which of should be lowered. */
/* NCFNEW : Degree +1 of the resulting polynom. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* YCVMAX : Auxiliary Table (max error on each dimension).
*/
/* ERRMAX : Precision of the approximation. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* > */
/* ***********************************************************************
*/
/* ---------------- Table of maximums of ((1-t2)3)*Ji(t) ---------------
*/
/* Parameter adjustments */
--ycvmax;
crvjac_dim1 = *ncofmx;
crvjac_offset = crvjac_dim1 + 1;
crvjac -= crvjac_offset;
/* Function Body */
/* ------------------- Init for error calculation -----------------------
*/
i__1 = *ndimen;
for (ii = 1; ii <= i__1; ++ii) {
ycvmax[ii] = 0.;
/* L100: */
}
/* ------ Min Degree that can be attained : Stop at 3 or NCFNEW ------
*/
idec = 7;
/* Computing MAX */
i__1 = idec, i__2 = *ncfnew + 1;
ncut = advapp_max(i__1,i__2);
/* -------------- Removal of high degree coefficients -----------
*/
/* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
*/
i__1 = *ncoeff;
for (ii = ncut; ii <= i__1; ++ii) {
/* Factor of renormalization. */
bidon = xmaxj[ii - idec];
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
bidon;
/* L310: */
}
/* L300: */
}
/* -------------- The error is the norm of the vector error ---------------
*/
*errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
/* --------------------------------- END --------------------------------
*/
return 0;
} /* mmaper6_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmaperx_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmaperx_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
integer *iordre,
doublereal *crvjac,
integer *ncfnew,
doublereal *ycvmax,
doublereal *errmax,
integer *iercod)
{
/* System generated locals */
integer crvjac_dim1, crvjac_offset;
/* Local variables */
integer jord;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Calculate the max. error of approximation made when */
/* only first NCFNEW coefficients of a curve are preserved
*/
/* degree NCOEFF-1 is written in the base of Jacobi of order IORDRE. */
/* KEYWORDS : */
/* ----------- */
/* JACOBI,LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX : Max. degree of the curve. */
/* NDIMEN : Space dimension. */
/* NCOEFF : Degree +1 of the curve. */
/* IORDRE : Order of continuity at the extremities. */
/* CRVJAC : Curve the degree which of should be lowered. */
/* NCFNEW : Degree +1 of the resulting polynom. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* YCVMAX : Auxiliary Table (max error on each dimension).
*/
/* ERRMAX : Precision of the approximation. */
/* IERCOD = 0, OK */
/* = 1, order of constraints (IORDRE) is not within the */
/* autorized values. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* Canceled and replaced MMAPERR. */
/* ***********************************************************************
*/
/* Parameter adjustments */
--ycvmax;
crvjac_dim1 = *ncofmx;
crvjac_offset = crvjac_dim1 + 1;
crvjac -= crvjac_offset;
/* Function Body */
*iercod = 0;
/* --> Order of Jacobi polynoms */
jord = ( *iordre + 1) << 1;
if (jord == 0) {
mmaper0_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
ycvmax[1], errmax);
} else if (jord == 2) {
mmaper2_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
ycvmax[1], errmax);
} else if (jord == 4) {
mmaper4_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
ycvmax[1], errmax);
} else if (jord == 6) {
mmaper6_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
ycvmax[1], errmax);
} else {
*iercod = 1;
}
/* ----------------------------------- Fin ------------------------------
*/
return 0;
} /* mmaperx_ */
//=======================================================================
//function : mmarc41_
//purpose :
//=======================================================================
int mmarc41_(integer *ndimax,
integer *ndimen,
integer *ncoeff,
doublereal *crvold,
doublereal *upara0,
doublereal *upara1,
doublereal *crvnew,
integer *iercod)
{
/* System generated locals */
integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
i__2, i__3;
/* Local variables */
integer nboct;
doublereal tbaux[61];
integer nd;
doublereal bid;
integer ncf, ncj;
/* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
/* IMPLICIT INTEGER (I-N) */
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Creation of curve C2(v) defined on (0,1) identic to */
/* curve C1(u) defined on (U0,U1) (change of parameter */
/* of a curve). */
/* KEYWORDS : */
/* ----------- */
/* LIMITATION, RESTRICTION, CURVE */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMAX : Space Dimensioning. */
/* NDIMEN : Curve Dimension. */
/* NCOEFF : Nb of coefficients of the curve. */
/* CRVOLD : Curve to be limited. */
/* UPARA0 : Min limit of the interval limiting the curve.
*/
/* UPARA1 : Max limit of the interval limiting the curve.
*/
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* CRVNEW : Relimited curve, defined on (0,1) and equal to */
/* CRVOLD defined on (U0,U1). */
/* IERCOD : = 0, OK */
/* =10, Nb of coeff. <1 or > 61. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Type Name */
/* MAERMSG MCRFILL MVCVIN2 */
/* MVCVINV */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ---> Algorithm used in this general case is based on the */
/* following principle : */
/* Let S(t) = a0 + a1*t + a2*t**2 + ... of degree NCOEFF-1, and */
/* U(t) = b0 + b1*t, then the coeff. of */
/* S(U(t)) are calculated step by step with help of table TBAUX. */
/* At each step number N (N=2 to NCOEFF), TBAUX(n) contains */
/* the n-th coefficient of U(t)**N for n=1 to N. (RBD) */
/* ---> Reference : KNUTH, 'The Art of Computer Programming', */
/* Vol. 2/'Seminumerical Algorithms', */
/* Ex. 11 p:451 et solution p:562. (RBD) */
/* ---> Removal of the input argument CRVOLD by CRVNEW is */
/* possible, which means that the call : */
/* CALL MMARC41(NDIMAX,NDIMEN,NCOEFF,CURVE,UPARA0,UPARA1 */
/* ,CURVE,IERCOD) */
/* is absolutely LEGAL. (RBD) */
/* > */
/* **********************************************************************
*/
/* Name of the routine */
/* Auxiliary table of coefficients of (UPARA1-UPARA0)T+UPARA0 */
/* with power N=1 to NCOEFF-1. */
/* Parameter adjustments */
crvnew_dim1 = *ndimax;
crvnew_offset = crvnew_dim1 + 1;
crvnew -= crvnew_offset;
crvold_dim1 = *ndimax;
crvold_offset = crvold_dim1 + 1;
crvold -= crvold_offset;
/* Function Body */
*iercod = 0;
/* **********************************************************************
*/
/* CASE WHEN PROCESSING CAN'T BE DONE */
/* **********************************************************************
*/
if (*ncoeff > 61 || *ncoeff < 1) {
*iercod = 10;
goto L9999;
}
/* **********************************************************************
*/
/* IF NO CHANGES */
/* **********************************************************************
*/
if (*ndimen == *ndimax && *upara0 == 0. && *upara1 == 1.) {
nboct = (*ndimax << 3) * *ncoeff;
AdvApp2Var_SysBase::mcrfill_(&nboct,
&crvold[crvold_offset],
&crvnew[crvnew_offset]);
goto L9999;
}
/* **********************************************************************
*/
/* INVERSION 3D : FAST PROCESSING */
/* **********************************************************************
*/
if (*upara0 == 1. && *upara1 == 0.) {
if (*ndimen == 3 && *ndimax == 3 && *ncoeff <= 21) {
mvcvinv_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
iercod);
goto L9999;
}
/* ******************************************************************
**** */
/* INVERSION 2D : FAST PROCESSING */
/* ******************************************************************
**** */
if (*ndimen == 2 && *ndimax == 2 && *ncoeff <= 21) {
mvcvin2_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
iercod);
goto L9999;
}
}
/* **********************************************************************
*/
/* GENERAL PROCESSING */
/* **********************************************************************
*/
/* -------------------------- Initializations ---------------------------
*/
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
crvnew[nd + crvnew_dim1] = crvold[nd + crvold_dim1];
/* L100: */
}
if (*ncoeff == 1) {
goto L9999;
}
tbaux[0] = *upara0;
tbaux[1] = *upara1 - *upara0;
/* ----------------------- Calculation of coeff. of CRVNEW ------------------
*/
i__1 = *ncoeff - 1;
for (ncf = 2; ncf <= i__1; ++ncf) {
/* ------------ Take into account NCF-th coeff. of CRVOLD --------
---- */
i__2 = ncf - 1;
for (ncj = 1; ncj <= i__2; ++ncj) {
bid = tbaux[ncj - 1];
i__3 = *ndimen;
for (nd = 1; nd <= i__3; ++nd) {
crvnew[nd + ncj * crvnew_dim1] += crvold[nd + ncf *
crvold_dim1] * bid;
/* L400: */
}
/* L300: */
}
bid = tbaux[ncf - 1];
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
crvnew[nd + ncf * crvnew_dim1] = crvold[nd + ncf * crvold_dim1] *
bid;
/* L500: */
}
/* --------- Calculate (NCF+1) coeff. of ((U1-U0)*t + U0)**(NCF) ---
---- */
bid = *upara1 - *upara0;
tbaux[ncf] = tbaux[ncf - 1] * bid;
for (ncj = ncf; ncj >= 2; --ncj) {
tbaux[ncj - 1] = tbaux[ncj - 1] * *upara0 + tbaux[ncj - 2] * bid;
/* L600: */
}
tbaux[0] *= *upara0;
/* L200: */
}
/* -------------- Take into account the last coeff. of CRVOLD -----------
*/
i__1 = *ncoeff - 1;
for (ncj = 1; ncj <= i__1; ++ncj) {
bid = tbaux[ncj - 1];
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
crvnew[nd + ncj * crvnew_dim1] += crvold[nd + *ncoeff *
crvold_dim1] * bid;
/* L800: */
}
/* L700: */
}
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
crvnew[nd + *ncoeff * crvnew_dim1] = crvold[nd + *ncoeff *
crvold_dim1] * tbaux[*ncoeff - 1];
/* L900: */
}
/* ---------------------------- The end ---------------------------------
*/
L9999:
if (*iercod != 0) {
AdvApp2Var_SysBase::maermsg_("MMARC41", iercod, 7L);
}
return 0 ;
} /* mmarc41_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmarcin_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmarcin_(integer *ndimax,
integer *ndim,
integer *ncoeff,
doublereal *crvold,
doublereal *u0,
doublereal *u1,
doublereal *crvnew,
integer *iercod)
{
/* System generated locals */
integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
i__2, i__3;
doublereal d__1;
/* Local variables */
doublereal x0, x1;
integer nd;
doublereal tabaux[61];
integer ibb;
doublereal bid;
integer ncf;
integer ncj;
doublereal eps3;
/* **********************************************************************
*//* FUNCTION : */
/* ---------- */
/* Creation of curve C2(v) defined on [U0,U1] identic to */
/* curve C1(u) defined on [-1,1] (change of parameter */
/* of a curve) with INVERSION of indices of the resulting table. */
/* KEYWORDS : */
/* ----------- */
/* GENERALIZED LIMITATION, RESTRICTION, INVERSION, CURVE */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMAX : Maximum Space Dimensioning. */
/* NDIMEN : Curve Dimension. */
/* NCOEFF : Nb of coefficients of the curve. */
/* CRVOLD : Curve to be limited. */
/* U0 : Min limit of the interval limiting the curve.
*/
/* U1 : Max limit of the interval limiting the curve.
*/
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* CRVNEW : Relimited curve, defined on [U0,U1] and equal to */
/* CRVOLD defined on [-1,1]. */
/* IERCOD : = 0, OK */
/* =10, Nb of coeff. <1 or > 61. */
/* =13, the requested interval of variation is null. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ---------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* **********************************************************************
*/
/* Name of the routine */
/* Auxiliary table of coefficients of X1*T+X0 */
/* with power N=1 to NCOEFF-1. */
/* Parameter adjustments */
crvnew_dim1 = *ndimax;
crvnew_offset = crvnew_dim1 + 1;
crvnew -= crvnew_offset;
crvold_dim1 = *ncoeff;
crvold_offset = crvold_dim1 + 1;
crvold -= crvold_offset;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 2) {
AdvApp2Var_SysBase::mgenmsg_("MMARCIN", 7L);
}
/* At zero machine it is tested if the output interval is not null */
AdvApp2Var_MathBase::mmveps3_(&eps3);
if ((d__1 = *u1 - *u0, advapp_abs(d__1)) < eps3) {
*iercod = 13;
goto L9999;
}
*iercod = 0;
/* **********************************************************************
*/
/* CASE WHEN THE PROCESSING IS IMPOSSIBLE */
/* **********************************************************************
*/
if (*ncoeff > 61 || *ncoeff < 1) {
*iercod = 10;
goto L9999;
}
/* **********************************************************************
*/
/* IF NO CHANGE OF THE INTERVAL OF DEFINITION */
/* (ONLY INVERSION OF INDICES OF TABLE CRVOLD) */
/* **********************************************************************
*/
if (*ndim == *ndimax && *u0 == -1. && *u1 == 1.) {
AdvApp2Var_MathBase::mmcvinv_(ndim, ncoeff, ndim, &crvold[crvold_offset], &crvnew[
crvnew_offset]);
goto L9999;
}
/* **********************************************************************
*/
/* CASE WHEN THE NEW INTERVAL OF DEFINITION IS [0,1] */
/* **********************************************************************
*/
if (*u0 == 0. && *u1 == 1.) {
mmcvstd_(ncoeff, ndimax, ncoeff, ndim, &crvold[crvold_offset], &
crvnew[crvnew_offset]);
goto L9999;
}
/* **********************************************************************
*/
/* GENERAL PROCESSING */
/* **********************************************************************
*/
/* -------------------------- Initialization ---------------------------
*/
x0 = -(*u1 + *u0) / (*u1 - *u0);
x1 = 2. / (*u1 - *u0);
i__1 = *ndim;
for (nd = 1; nd <= i__1; ++nd) {
crvnew[nd + crvnew_dim1] = crvold[nd * crvold_dim1 + 1];
/* L100: */
}
if (*ncoeff == 1) {
goto L9999;
}
tabaux[0] = x0;
tabaux[1] = x1;
/* ----------------------- Calculation of coeff. of CRVNEW ------------------
*/
i__1 = *ncoeff - 1;
for (ncf = 2; ncf <= i__1; ++ncf) {
/* ------------ Take into account the NCF-th coeff. of CRVOLD --------
---- */
i__2 = ncf - 1;
for (ncj = 1; ncj <= i__2; ++ncj) {
bid = tabaux[ncj - 1];
i__3 = *ndim;
for (nd = 1; nd <= i__3; ++nd) {
crvnew[nd + ncj * crvnew_dim1] += crvold[ncf + nd *
crvold_dim1] * bid;
/* L400: */
}
/* L300: */
}
bid = tabaux[ncf - 1];
i__2 = *ndim;
for (nd = 1; nd <= i__2; ++nd) {
crvnew[nd + ncf * crvnew_dim1] = crvold[ncf + nd * crvold_dim1] *
bid;
/* L500: */
}
/* --------- Calculation of (NCF+1) coeff. of [X1*t + X0]**(NCF) --------
---- */
tabaux[ncf] = tabaux[ncf - 1] * x1;
for (ncj = ncf; ncj >= 2; --ncj) {
tabaux[ncj - 1] = tabaux[ncj - 1] * x0 + tabaux[ncj - 2] * x1;
/* L600: */
}
tabaux[0] *= x0;
/* L200: */
}
/* -------------- Take into account the last coeff. of CRVOLD -----------
*/
i__1 = *ncoeff - 1;
for (ncj = 1; ncj <= i__1; ++ncj) {
bid = tabaux[ncj - 1];
i__2 = *ndim;
for (nd = 1; nd <= i__2; ++nd) {
crvnew[nd + ncj * crvnew_dim1] += crvold[*ncoeff + nd *
crvold_dim1] * bid;
/* L800: */
}
/* L700: */
}
i__1 = *ndim;
for (nd = 1; nd <= i__1; ++nd) {
crvnew[nd + *ncoeff * crvnew_dim1] = crvold[*ncoeff + nd *
crvold_dim1] * tabaux[*ncoeff - 1];
/* L900: */
}
/* ---------------------------- The end ---------------------------------
*/
L9999:
if (*iercod > 0) {
AdvApp2Var_SysBase::maermsg_("MMARCIN", iercod, 7L);
}
if (ibb >= 2) {
AdvApp2Var_SysBase::mgsomsg_("MMARCIN", 7L);
}
return 0;
} /* mmarcin_ */
//=======================================================================
//function : mmatvec_
//purpose :
//=======================================================================
int mmatvec_(integer *nligne,
integer *,//ncolon,
integer *gposit,
integer *,//gnstoc,
doublereal *gmatri,
doublereal *vecin,
integer *deblig,
doublereal *vecout,
integer *iercod)
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
logical ldbg;
integer jmin, jmax, i__, j, k;
doublereal somme;
integer aux;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Produce vector matrix in form of profile */
/* MOTS CLES : */
/* ----------- */
/* RESERVE, MATRIX, PRODUCT, VECTOR, PROFILE */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* NLIGNE : Line number of the matrix of constraints */
/* NCOLON : Number of column of the matrix of constraints */
/* GNSTOC: Number of coefficients in the profile of matrix GMATRI */
/* GPOSIT: Table of positioning of terms of storage */
/* GPOSIT(1,I) contains the number of terms-1 on the line I */
/* in the profile of the matrix. */
/* GPOSIT(2,I) contains the index of storage of diagonal term*/
/* of line I */
/* GPOSIT(3,I) contains the index of column of the first term of */
/* profile of line I */
/* GNSTOC: Number of coefficients in the profile of matrix */
/* GMATRI */
/* GMATRI : Matrix of constraints in form of profile */
/* VECIN : Input vector */
/* DEBLIG : Line indexusing which the vector matrix is calculated */
/**/
/* OUTPUT ARGUMENTS */
/* --------------------- */
/* VECOUT : VECTOR PRODUCT */
/* IERCOD : ERROR CODE */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* INITIALISATIONS */
/* ***********************************************************************
*/
/* Parameter adjustments */
--vecout;
gposit -= 4;
--vecin;
--gmatri;
/* Function Body */
ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
if (ldbg) {
AdvApp2Var_SysBase::mgenmsg_("MMATVEC", 7L);
}
*iercod = 0;
/* ***********************************************************************
*/
/* Processing */
/* ***********************************************************************
*/
AdvApp2Var_SysBase::mvriraz_(nligne,
&vecout[1]);
i__1 = *nligne;
for (i__ = *deblig; i__ <= i__1; ++i__) {
somme = 0.;
jmin = gposit[i__ * 3 + 3];
jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
i__2 = jmax;
for (j = jmin; j <= i__2; ++j) {
k = j + aux;
somme += gmatri[k] * vecin[j];
}
vecout[i__] = somme;
}
goto L9999;
/* ***********************************************************************
*/
/* ERROR PROCESSING */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
L9999:
/* ___ DESALLOCATION, ... */
AdvApp2Var_SysBase::maermsg_("MMATVEC", iercod, 7L);
if (ldbg) {
AdvApp2Var_SysBase::mgsomsg_("MMATVEC", 7L);
}
return 0 ;
} /* mmatvec_ */
//=======================================================================
//function : mmbulld_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmbulld_(integer *nbcoln,
integer *nblign,
doublereal *dtabtr,
integer *numcle)
{
/* System generated locals */
integer dtabtr_dim1, dtabtr_offset, i__1, i__2;
/* Local variables */
logical ldbg;
doublereal daux;
integer nite1, nite2, nchan, i1, i2;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Parsing of columns of a table of integers in increasing order */
/* KEYWORDS : */
/* ----------- */
/* POINT-ENTRY, PARSING */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* - NBCOLN : Number of columns in the table */
/* - NBLIGN : Number of lines in the table */
/* - DTABTR : Table of integers to be parsed */
/* - NUMCLE : Position of the key on the column */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* - DTABTR : Parsed table */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* Particularly performant if the table is almost parsed */
/* In the opposite case it is better to use MVSHELD */
/* ***********************************************************************
*/
/* Parameter adjustments */
dtabtr_dim1 = *nblign;
dtabtr_offset = dtabtr_dim1 + 1;
dtabtr -= dtabtr_offset;
/* Function Body */
ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
if (ldbg) {
AdvApp2Var_SysBase::mgenmsg_("MMBULLD", 7L);
}
nchan = 1;
nite1 = *nbcoln;
nite2 = 2;
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
/* ---->ALGORITHM in N^2 / 2 additional iteration */
while(nchan != 0) {
/* ----> Parsing from left to the right */
nchan = 0;
i__1 = nite1;
for (i1 = nite2; i1 <= i__1; ++i1) {
if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1 - 1)
* dtabtr_dim1]) {
i__2 = *nblign;
for (i2 = 1; i2 <= i__2; ++i2) {
daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
dtabtr_dim1];
dtabtr[i2 + i1 * dtabtr_dim1] = daux;
}
if (nchan == 0) {
nchan = 1;
}
}
}
--nite1;
/* ----> Parsing from right to the left */
if (nchan != 0) {
nchan = 0;
i__1 = nite2;
for (i1 = nite1; i1 >= i__1; --i1) {
if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1
- 1) * dtabtr_dim1]) {
i__2 = *nblign;
for (i2 = 1; i2 <= i__2; ++i2) {
daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
dtabtr_dim1];
dtabtr[i2 + i1 * dtabtr_dim1] = daux;
}
if (nchan == 0) {
nchan = 1;
}
}
}
++nite2;
}
}
goto L9999;
/* ***********************************************************************
*/
/* ERROR PROCESSING */
/* ***********************************************************************
*/
/* ----> No errors at calling functions, only tests and loops. */
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
L9999:
if (ldbg) {
AdvApp2Var_SysBase::mgsomsg_("MMBULLD", 7L);
}
return 0 ;
} /* mmbulld_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmcdriv_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmcdriv_(integer *ndimen,
integer *ncoeff,
doublereal *courbe,
integer *ideriv,
integer *ncofdv,
doublereal *crvdrv)
{
/* System generated locals */
integer courbe_dim1, courbe_offset, crvdrv_dim1, crvdrv_offset, i__1,
i__2;
/* Local variables */
integer i__, j, k;
doublereal mfactk, bid;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Calculate matrix of a derivate curve of order IDERIV. */
/* with input parameters other than output parameters. */
/* KEYWORDS : */
/* ----------- */
/* COEFFICIENTS,CURVE,DERIVATE I-EME. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMEN : Space dimension (2 or 3 in general) */
/* NCOEFF : Degree +1 of the curve. */
/* COURBE : Table of coefficients of the curve. */
/* IDERIV : Required order of derivation : 1=1st derivate, etc... */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* NCOFDV : Degree +1 of the derivative of order IDERIV of the curve. */
/* CRVDRV : Table of coefficients of the derivative of order IDERIV */
/* of the curve. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ---> It is possible to take as output argument the curve */
/* and the number of coeff passed at input by making : */
/* CALL MMCDRIV(NDIMEN,NCOEFF,COURBE,IDERIV,NCOEFF,COURBE). */
/* After this call, NCOEFF does the number of coeff of the derived */
/* curve the coefficients which of are stored in CURVE. */
/* Attention to the coefficients of CURVE of rank superior to */
/* NCOEFF : they are not set to zero. */
/* ---> Algorithm : */
/* The code below was written basing on the following algorithm:
*/
/* Let P(t) = a1 + a2*t + ... an*t**n. Derivate of order k of P */
/* (containing n-k coefficients) is calculated as follows : */
/* Pk(t) = a(k+1)*CNP(k,k)*k! */
/* + a(k+2)*CNP(k+1,k)*k! * t */
/* . */
/* . */
/* . */
/* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
/* ***********************************************************************
*/
/* -------------- Case when the order of derivative is -------------------
*/
/* ---------------- greater than the degree of the curve ---------------------
*/
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Serves to provide the coefficients of binome (Pascal's triangle). */
/* KEYWORDS : */
/* ----------- */
/* Binomial coeff from 0 to 60. read only . init par block data */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* Binomial coefficients form a triangular matrix. */
/* This matrix is completed in table CNP by its transposition. */
/* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
/* Initialization is done by block-data MMLLL09.RES, */
/* created by program MQINICNP.FOR). */
/* **********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
crvdrv_dim1 = *ndimen;
crvdrv_offset = crvdrv_dim1 + 1;
crvdrv -= crvdrv_offset;
courbe_dim1 = *ndimen;
courbe_offset = courbe_dim1 + 1;
courbe -= courbe_offset;
/* Function Body */
if (*ideriv >= *ncoeff) {
i__1 = *ndimen;
for (i__ = 1; i__ <= i__1; ++i__) {
crvdrv[i__ + crvdrv_dim1] = 0.;
/* L10: */
}
*ncofdv = 1;
goto L9999;
}
/* **********************************************************************
*/
/* General processing */
/* **********************************************************************
*/
/* --------------------- Calculation of Factorial(IDERIV) ------------------
*/
k = *ideriv;
mfactk = 1.;
i__1 = k;
for (i__ = 2; i__ <= i__1; ++i__) {
mfactk *= i__;
/* L50: */
}
/* ------------ Calculation of coeff of the derived of order IDERIV ----------
*/
/* ---> Attention : coefficient binomial C(n,m) is represented in */
/* MCCNP by CNP(N+1,M+1). */
i__1 = *ncoeff;
for (j = k + 1; j <= i__1; ++j) {
bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
i__2 = *ndimen;
for (i__ = 1; i__ <= i__2; ++i__) {
crvdrv[i__ + (j - k) * crvdrv_dim1] = bid * courbe[i__ + j *
courbe_dim1];
/* L200: */
}
/* L100: */
}
*ncofdv = *ncoeff - *ideriv;
/* -------------------------------- The end -----------------------------
*/
L9999:
return 0;
} /* mmcdriv_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmcglc1_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmcglc1_(integer *ndimax,
integer *ndimen,
integer *ncoeff,
doublereal *courbe,
doublereal *tdebut,
doublereal *tfinal,
doublereal *epsiln,
doublereal *xlongc,
doublereal *erreur,
integer *iercod)
{
/* System generated locals */
integer courbe_dim1, courbe_offset, i__1;
doublereal d__1;
/* Local variables */
integer ndec;
doublereal tdeb, tfin;
integer iter;
doublereal oldso = 0.;
integer itmax;
doublereal sottc;
integer kk, ibb;
doublereal dif, pas;
doublereal som;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Allows calculating the length of an arc of curve POLYNOMIAL */
/* on an interval [A,B]. */
/* KEYWORDS : */
/* ----------- */
/* LENGTH,CURVE,GAUSS,PRIVATE. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMAX : Max. number of lines of tables */
/* (i.e. max. nb of polynoms). */
/* NDIMEN : Dimension of the space (nb of polynoms). */
/* NCOEFF : Nb of coefficients of the polynom. This is degree + 1.
*/
/* COURBE(NDIMAX,NCOEFF) : Coefficients of the curve. */
/* TDEBUT : Lower limit of the interval of integration for */
/* length calculation. */
/* TFINAL : Upper limit of the interval of integration for */
/* length calculation. */
/* EPSILN : REQIRED precision for length calculation. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* XLONGC : Length of the arc of curve */
/* ERREUR : Precision OBTAINED for the length calculation. */
/* IERCOD : Error code, 0 OK, >0 Serious error. */
/* = 1 Too much iterations, the best calculated resultat */
/* (is almost ERROR) */
/* = 2 Pb MMLONCV (no result) */
/* = 3 NDIM or NCOEFF invalid (no result) */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* The polynom is actually a set of polynoms with */
/* coefficients arranged in a table of 2 indices, */
/* each line relative to the polynom. */
/* The polynom is defined by these coefficients ordered */
/* by increasing power of the variable. */
/* All polynoms have the same number of coefficients (the */
/* same degree). */
/* This program cancels and replaces LENGCV, MLONGC and MLENCV. */
/* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
/* > */
/* ***********************************************************************
*/
/* Name of the routine */
/* ------------------------ General Initialization ---------------------
*/
/* Parameter adjustments */
courbe_dim1 = *ndimax;
courbe_offset = courbe_dim1 + 1;
courbe -= courbe_offset;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 2) {
AdvApp2Var_SysBase::mgenmsg_("MMCGLC1", 7L);
}
*iercod = 0;
*xlongc = 0.;
*erreur = 0.;
/* ------ Test of equity of limits */
if (*tdebut == *tfinal) {
*iercod = 0;
goto L9999;
}
/* ------ Test of the dimension and the number of coefficients */
if (*ndimen <= 0 || *ncoeff <= 0) {
goto L9003;
}
/* ----- Nb of current cutting, nb of iteration, */
/* max nb of iterations */
ndec = 1;
iter = 1;
itmax = 13;
/* ------ Variation of the nb of intervals */
/* Multiplied by 2 at each iteration */
L5000:
pas = (*tfinal - *tdebut) / ndec;
sottc = 0.;
/* ------ Loop on all current NDEC intervals */
i__1 = ndec;
for (kk = 1; kk <= i__1; ++kk) {
/* ------ Limits of the current integration interval */
tdeb = *tdebut + (kk - 1) * pas;
tfin = tdeb + pas;
mmloncv_(ndimax, ndimen, ncoeff, &courbe[courbe_offset], &tdeb, &tfin,
&som, iercod);
if (*iercod > 0) {
goto L9002;
}
sottc += som;
/* L100: */
}
/* ----------------- Test of the maximum number of iterations ------------
*/
/* Test if passes at least once ** */
if (iter == 1) {
oldso = sottc;
ndec <<= 1;
++iter;
goto L5000;
} else {
/* ------ Take into account DIF - Test of convergence */
++iter;
dif = (d__1 = sottc - oldso, advapp_abs(d__1));
/* ------ If DIF is OK, leave..., otherwise: */
if (dif > *epsiln) {
/* ------ If nb iteration exceeded, leave */
if (iter > itmax) {
*iercod = 1;
goto L9000;
} else {
/* ------ Otherwise continue by cutting the initial interval.
*/
oldso = sottc;
ndec <<= 1;
goto L5000;
}
}
}
/* ------------------------------ THE END -------------------------------
*/
L9000:
*xlongc = sottc;
*erreur = dif;
goto L9999;
/* ---> PB in MMLONCV */
L9002:
*iercod = 2;
goto L9999;
/* ---> NCOEFF or NDIM invalid. */
L9003:
*iercod = 3;
goto L9999;
L9999:
if (*iercod > 0) {
AdvApp2Var_SysBase::maermsg_("MMCGLC1", iercod, 7L);
}
if (ibb >= 2) {
AdvApp2Var_SysBase::mgsomsg_("MMCGLC1", 7L);
}
return 0;
} /* mmcglc1_ */
//=======================================================================
//function : mmchole_
//purpose :
//=======================================================================
int mmchole_(integer *,//mxcoef,
integer *dimens,
doublereal *amatri,
integer *aposit,
integer *posuiv,
doublereal *chomat,
integer *iercod)
{
/* System generated locals */
integer i__1, i__2, i__3;
doublereal d__1;
/* Builtin functions */
//double sqrt();
/* Local variables */
logical ldbg;
integer kmin, i__, j, k;
doublereal somme;
integer ptini, ptcou;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- T */
/* Produce decomposition of choleski of matrix A in S.S */
/* Calculate inferior triangular matrix S. */
/* KEYWORDS : */
/* ----------- */
/* RESOLUTION, MFACTORISATION, MATRIX_PROFILE, CHOLESKI */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* MXCOEF : Max number of terms in the hessian profile */
/* DIMENS : Dimension of the problem */
/* AMATRI(MXCOEF) : Coefficients of the matrix profile */
/* APOSIT(1,*) : Distance diagonal-left extremity of the line
*/
/* APOSIT(2,*) : Position of diagonal terms in HESSIE */
/* POSUIV(MXCOEF) : first line inferior not out of profile */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* CHOMAT(MXCOEF) : Inferior triangular matrix preserving the */
/* profile of AMATRI. */
/* IERCOD : error code */
/* = 0 : ok */
/* = 1 : non-defined positive matrix */
/* COMMONS USED : */
/* ------------------ */
/* .Neant. */
/* REFERENCES CALLED : */
/* ---------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* DEBUG LEVEL = 4 */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* INITIALISATIONS */
/* ***********************************************************************
*/
/* Parameter adjustments */
--chomat;
--posuiv;
--amatri;
aposit -= 3;
/* Function Body */
ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
if (ldbg) {
AdvApp2Var_SysBase::mgenmsg_("MMCHOLE", 7L);
}
*iercod = 0;
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
i__1 = *dimens;
for (j = 1; j <= i__1; ++j) {
ptini = aposit[(j << 1) + 2];
somme = 0.;
i__2 = ptini - 1;
for (k = ptini - aposit[(j << 1) + 1]; k <= i__2; ++k) {
/* Computing 2nd power */
d__1 = chomat[k];
somme += d__1 * d__1;
}
if (amatri[ptini] - somme < 1e-32) {
goto L9101;
}
chomat[ptini] = sqrt(amatri[ptini] - somme);
ptcou = ptini;
while(posuiv[ptcou] > 0) {
i__ = posuiv[ptcou];
ptcou = aposit[(i__ << 1) + 2] - (i__ - j);
/* Calculate the sum of S .S for k =1 a j-1 */
/* ik jk */
somme = 0.;
/* Computing MAX */
i__2 = i__ - aposit[(i__ << 1) + 1], i__3 = j - aposit[(j << 1) +
1];
kmin = advapp_max(i__2,i__3);
i__2 = j - 1;
for (k = kmin; k <= i__2; ++k) {
somme += chomat[aposit[(i__ << 1) + 2] - (i__ - k)] * chomat[
aposit[(j << 1) + 2] - (j - k)];
}
chomat[ptcou] = (amatri[ptcou] - somme) / chomat[ptini];
}
}
goto L9999;
/* ***********************************************************************
*/
/* ERROR PROCESSING */
/* ***********************************************************************
*/
L9101:
*iercod = 1;
goto L9999;
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
L9999:
AdvApp2Var_SysBase::maermsg_("MMCHOLE", iercod, 7L);
if (ldbg) {
AdvApp2Var_SysBase::mgsomsg_("MMCHOLE", 7L);
}
return 0 ;
} /* mmchole_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmcvctx_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmcvctx_(integer *ndimen,
integer *ncofmx,
integer *nderiv,
doublereal *ctrtes,
doublereal *crvres,
doublereal *tabaux,
doublereal *xmatri,
integer *iercod)
{
/* System generated locals */
integer ctrtes_dim1, ctrtes_offset, crvres_dim1, crvres_offset,
xmatri_dim1, xmatri_offset, tabaux_dim1, tabaux_offset, i__1,
i__2;
/* Local variables */
integer moup1, nordr;
integer nd;
integer ibb, ncf, ndv;
doublereal eps1;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Calculate a polynomial curve checking the */
/* passage constraints (interpolation) */
/* from first derivatives, etc... to extremities. */
/* Parameters at the extremities are supposed to be -1 and 1. */
/* KEYWORDS : */
/* ----------- */
/* ALL, AB_SPECIFI::CONSTRAINTS&,INTERPOLATION,&CURVE */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMEN : Space Dimension. */
/* NCOFMX : Nb of coeff. of curve CRVRES on each */
/* dimension. */
/* NDERIV : Order of constraint with derivatives : */
/* 0 --> interpolation simple. */
/* 1 --> interpolation+constraints with 1st. */
/* 2 --> cas (0)+ (1) + " " 2nd derivatives. */
/* etc... */
/* CTRTES : Table of constraints. */
/* CTRTES(*,1,*) = contraints at -1. */
/* CTRTES(*,2,*) = contraints at 1. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* CRVRES : Resulting curve defined on (-1,1). */
/* TABAUX : Auxilliary matrix. */
/* XMATRI : Auxilliary matrix. */
/* COMMONS UTILISES : */
/* ---------------- */
/* .Neant. */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Type Name */
/* MAERMSG R*8 DFLOAT MGENMSG */
/* MGSOMSG MMEPS1 MMRSLW */
/* I*4 MNFNDEB */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* The polynom (or the curve) is calculated by solving a */
/* system of linear equations. If the imposed degree is great */
/* it is preferable to call a routine based on */
/* Lagrange or Hermite interpolation depending on the case. */
/* (for a high degree the matrix of the system can be badly */
/* conditionned). */
/* This routine returns a curve defined in (-1,1). */
/* In general case, it is necessary to use MCVCTG. */
/* > */
/* ***********************************************************************
*/
/* Name of the routine */
/* Parameter adjustments */
crvres_dim1 = *ncofmx;
crvres_offset = crvres_dim1 + 1;
crvres -= crvres_offset;
xmatri_dim1 = *nderiv + 1;
xmatri_offset = xmatri_dim1 + 1;
xmatri -= xmatri_offset;
tabaux_dim1 = *nderiv + 1 + *ndimen;
tabaux_offset = tabaux_dim1 + 1;
tabaux -= tabaux_offset;
ctrtes_dim1 = *ndimen;
ctrtes_offset = ctrtes_dim1 * 3 + 1;
ctrtes -= ctrtes_offset;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 3) {
AdvApp2Var_SysBase::mgenmsg_("MMCVCTX", 7L);
}
/* Precision. */
AdvApp2Var_MathBase::mmeps1_(&eps1);
/* ****************** CALCULATION OF EVEN COEFFICIENTS *********************
*/
/* ------------------------- Initialization -----------------------------
*/
nordr = *nderiv + 1;
i__1 = nordr;
for (ncf = 1; ncf <= i__1; ++ncf) {
tabaux[ncf + tabaux_dim1] = 1.;
/* L100: */
}
/* ---------------- Calculation of terms corresponding to derivatives -------
*/
i__1 = nordr;
for (ndv = 2; ndv <= i__1; ++ndv) {
i__2 = nordr;
for (ncf = 1; ncf <= i__2; ++ncf) {
tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
tabaux_dim1] * (doublereal) ((ncf << 1) - ndv);
/* L300: */
}
/* L200: */
}
/* ------------------ Writing the second member -----------------------
*/
moup1 = 1;
i__1 = nordr;
for (ndv = 1; ndv <= i__1; ++ndv) {
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
+ 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
* ctrtes_dim1]) / 2.;
/* L500: */
}
moup1 = -moup1;
/* L400: */
}
/* -------------------- Resolution of the system ---------------------------
*/
mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
xmatri_offset], iercod);
if (*iercod > 0) {
goto L9999;
}
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
i__2 = nordr;
for (ncf = 1; ncf <= i__2; ++ncf) {
crvres[(ncf << 1) - 1 + nd * crvres_dim1] = xmatri[ncf + nd *
xmatri_dim1];
/* L700: */
}
/* L600: */
}
/* ***************** CALCULATION OF UNEVEN COEFFICIENTS ********************
*/
/* ------------------------- Initialization -----------------------------
*/
i__1 = nordr;
for (ncf = 1; ncf <= i__1; ++ncf) {
tabaux[ncf + tabaux_dim1] = 1.;
/* L1100: */
}
/* ---------------- Calculation of terms corresponding to derivatives -------
*/
i__1 = nordr;
for (ndv = 2; ndv <= i__1; ++ndv) {
i__2 = nordr;
for (ncf = 1; ncf <= i__2; ++ncf) {
tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
tabaux_dim1] * (doublereal) ((ncf << 1) - ndv + 1);
/* L1300: */
}
/* L1200: */
}
/* ------------------ Writing of the second member -----------------------
*/
moup1 = -1;
i__1 = nordr;
for (ndv = 1; ndv <= i__1; ++ndv) {
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
+ 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
* ctrtes_dim1]) / 2.;
/* L1500: */
}
moup1 = -moup1;
/* L1400: */
}
/* -------------------- Solution of the system ---------------------------
*/
mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
xmatri_offset], iercod);
if (*iercod > 0) {
goto L9999;
}
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
i__2 = nordr;
for (ncf = 1; ncf <= i__2; ++ncf) {
crvres[(ncf << 1) + nd * crvres_dim1] = xmatri[ncf + nd *
xmatri_dim1];
/* L1700: */
}
/* L1600: */
}
/* --------------------------- The end ----------------------------------
*/
L9999:
if (*iercod != 0) {
AdvApp2Var_SysBase::maermsg_("MMCVCTX", iercod, 7L);
}
if (ibb >= 3) {
AdvApp2Var_SysBase::mgsomsg_("MMCVCTX", 7L);
}
return 0 ;
} /* mmcvctx_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmcvinv_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmcvinv_(integer *ndimax,
integer *ncoef,
integer *ndim,
doublereal *curveo,
doublereal *curve)
{
/* Initialized data */
static char nomprg[8+1] = "MMCVINV ";
/* System generated locals */
integer curve_dim1, curve_offset, curveo_dim1, curveo_offset, i__1, i__2;
/* Local variables */
integer i__, nd, ibb;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Inversion of arguments of the final curve. */
/* KEYWORDS : */
/* ----------- */
/* SMOOTHING,CURVE */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIM: Space Dimension. */
/* NCOEF: Degree of the polynom. */
/* CURVEO: The curve before inversion. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* CURVE: The curve after inversion. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES APPELEES : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ***********************************************************************
*/
/* The name of the routine */
/* Parameter adjustments */
curve_dim1 = *ndimax;
curve_offset = curve_dim1 + 1;
curve -= curve_offset;
curveo_dim1 = *ncoef;
curveo_offset = curveo_dim1 + 1;
curveo -= curveo_offset;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 2) {
AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
}
i__1 = *ncoef;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = *ndim;
for (nd = 1; nd <= i__2; ++nd) {
curve[nd + i__ * curve_dim1] = curveo[i__ + nd * curveo_dim1];
/* L300: */
}
}
/* L9999: */
return 0;
} /* mmcvinv_ */
//=======================================================================
//function : mmcvstd_
//purpose :
//=======================================================================
int mmcvstd_(integer *ncofmx,
integer *ndimax,
integer *ncoeff,
integer *ndimen,
doublereal *crvcan,
doublereal *courbe)
{
/* System generated locals */
integer courbe_dim1, crvcan_dim1, crvcan_offset, i__1, i__2, i__3;
/* Local variables */
integer ndeg, i__, j, j1, nd, ibb;
doublereal bid;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Transform curve defined between [-1,1] into [0,1]. */
/* KEYWORDS : */
/* ----------- */
/* LIMITATION,RESTRICTION,CURVE */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMAX : Dimension of the space. */
/* NDIMEN : Dimension of the curve. */
/* NCOEFF : Degree of the curve. */
/* CRVCAN(NCOFMX,NDIMEN): The curve is defined at the interval [-1,1]. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* CURVE(NDIMAX,NCOEFF): Curve defined at the interval [0,1]. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* Name of the program. */
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Provides binomial coefficients (Pascal triangle). */
/* KEYWORDS : */
/* ----------- */
/* Binomial coefficient from 0 to 60. read only . init by block data */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* Binomial coefficients form a triangular matrix. */
/* This matrix is completed in table CNP by its transposition. */
/* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
/* Initialization is done with block-data MMLLL09.RES, */
/* created by the program MQINICNP.FOR. */
/* > */
/* **********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
courbe_dim1 = *ndimax;
--courbe;
crvcan_dim1 = *ncofmx;
crvcan_offset = crvcan_dim1;
crvcan -= crvcan_offset;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 3) {
AdvApp2Var_SysBase::mgenmsg_("MMCVSTD", 7L);
}
ndeg = *ncoeff - 1;
/* ------------------ Construction of the resulting curve ----------------
*/
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
i__2 = ndeg;
for (j = 0; j <= i__2; ++j) {
bid = 0.;
i__3 = ndeg;
for (i__ = j; i__ <= i__3; i__ += 2) {
bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
* 61];
/* L410: */
}
courbe[nd + j * courbe_dim1] = bid;
bid = 0.;
j1 = j + 1;
i__3 = ndeg;
for (i__ = j1; i__ <= i__3; i__ += 2) {
bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
* 61];
/* L420: */
}
courbe[nd + j * courbe_dim1] -= bid;
/* L400: */
}
/* L300: */
}
/* ------------------- Renormalization of the CURVE -------------------------
*/
bid = 1.;
i__1 = ndeg;
for (i__ = 0; i__ <= i__1; ++i__) {
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
courbe[nd + i__ * courbe_dim1] *= bid;
/* L510: */
}
bid *= 2.;
/* L500: */
}
/* ----------------------------- The end --------------------------------
*/
if (ibb >= 3) {
AdvApp2Var_SysBase::mgsomsg_("MMCVSTD", 7L);
}
return 0;
} /* mmcvstd_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmdrc11_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmdrc11_(integer *iordre,
integer *ndimen,
integer *ncoeff,
doublereal *courbe,
doublereal *points,
doublereal *mfactab)
{
/* System generated locals */
integer courbe_dim1, courbe_offset, points_dim2, points_offset, i__1,
i__2;
/* Local variables */
integer ndeg, i__, j, ndgcb, nd, ibb;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Calculation of successive derivatives of equation CURVE with */
/* parameters -1, 1 from order 0 to order IORDRE */
/* included. The calculation is produced without knowing the coefficients of */
/* derivatives of the curve. */
/* KEYWORDS : */
/* ----------- */
/* POSITIONING,EXTREMITIES,CURVE,DERIVATIVE. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* IORDRE : Maximum order of calculation of derivatives. */
/* NDIMEN : Dimension of the space. */
/* NCOEFF : Number of coefficients of the curve (degree+1). */
/* COURBE : Table of coefficients of the curve. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* POINTS : Table of values of consecutive derivatives */
/* of parameters -1.D0 and 1.D0. */
/* MFACTAB : Auxiliary table for calculation of factorial(I).
*/
/* COMMONS USED : */
/* ---------------- */
/* None. */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ---> ATTENTION, the coefficients of the curve are */
/* in a reverse order. */
/* ---> The algorithm of calculation of derivatives is based on */
/* generalization of Horner scheme : */
/* k 2 */
/* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
/* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
/* aj = a(j-1).x + u(k-j) */
/* bj = b(j-1).x + a(j-1) */
/* cj = c(j-1).x + b(j-1) */
/* So : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
/* The algorithm is generalized easily for calculation of */
/* (n) */
/* C (x) . */
/* --------- */
/* n! */
/* Reference : D. KNUTH, "The Art of Computer Programming" */
/* --------- Vol. 2/Seminumerical Algorithms */
/* Addison-Wesley Pub. Co. (1969) */
/* pages 423-425. */
/* > */
/* **********************************************************************
*/
/* Name of the routine */
/* Parameter adjustments */
points_dim2 = *iordre + 1;
points_offset = (points_dim2 << 1) + 1;
points -= points_offset;
courbe_dim1 = *ncoeff;
courbe_offset = courbe_dim1;
courbe -= courbe_offset;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 2) {
AdvApp2Var_SysBase::mgenmsg_("MMDRC11", 7L);
}
if (*iordre < 0 || *ncoeff < 1) {
goto L9999;
}
/* ------------------- Initialization of table POINTS -----------------
*/
ndgcb = *ncoeff - 1;
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
points[(nd * points_dim2 << 1) + 1] = courbe[ndgcb + nd * courbe_dim1]
;
points[(nd * points_dim2 << 1) + 2] = courbe[ndgcb + nd * courbe_dim1]
;
/* L100: */
}
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
i__2 = *iordre;
for (j = 1; j <= i__2; ++j) {
points[((j + nd * points_dim2) << 1) + 1] = 0.;
points[((j + nd * points_dim2) << 1) + 2] = 0.;
/* L400: */
}
/* L300: */
}
/* Calculation with parameter -1 and 1 */
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
i__2 = ndgcb;
for (ndeg = 1; ndeg <= i__2; ++ndeg) {
for (i__ = *iordre; i__ >= 1; --i__) {
points[((i__ + nd * points_dim2) << 1) + 1] = -points[((i__ + nd
* points_dim2) << 1) + 1] + points[((i__ - 1 + nd *
points_dim2) << 1) + 1];
points[((i__ + nd * points_dim2) << 1) + 2] += points[((i__ - 1
+ nd * points_dim2) << 1) + 2];
/* L800: */
}
points[(nd * points_dim2 << 1) + 1] = -points[(nd * points_dim2 <<
1) + 1] + courbe[ndgcb - ndeg + nd * courbe_dim1];
points[(nd * points_dim2 << 1) + 2] += courbe[ndgcb - ndeg + nd *
courbe_dim1];
/* L700: */
}
/* L600: */
}
/* --------------------- Multiplication by factorial(I) --------------
*/
if (*iordre > 1) {
mfac_(&mfactab[1], iordre);
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
i__2 = *iordre;
for (i__ = 2; i__ <= i__2; ++i__) {
points[((i__ + nd * points_dim2) << 1) + 1] = mfactab[i__] *
points[((i__ + nd * points_dim2) << 1) + 1];
points[((i__ + nd * points_dim2) << 1) + 2] = mfactab[i__] *
points[((i__ + nd * points_dim2) << 1) + 2];
/* L1000: */
}
/* L900: */
}
}
/* ---------------------------- End -------------------------------------
*/
L9999:
if (ibb >= 2) {
AdvApp2Var_SysBase::mgsomsg_("MMDRC11", 7L);
}
return 0;
} /* mmdrc11_ */
//=======================================================================
//function : mmdrvcb_
//purpose :
//=======================================================================
int mmdrvcb_(integer *ideriv,
integer *ndim,
integer *ncoeff,
doublereal *courbe,
doublereal *tparam,
doublereal *tabpnt,
integer *iercod)
{
/* System generated locals */
integer courbe_dim1, tabpnt_dim1, i__1, i__2, i__3;
/* Local variables */
integer ndeg, i__, j, nd, ndgcrb, iptpnt, ibb;
/* *********************************************************************** */
/* FUNCTION : */
/* ---------- */
/* Calculation of successive derivatives of equation CURVE with */
/* parameter TPARAM from order 0 to order IDERIV included. */
/* The calculation is produced without knowing the coefficients of */
/* derivatives of the CURVE. */
/* KEYWORDS : */
/* ----------- */
/* POSITIONING,PARAMETER,CURVE,DERIVATIVE. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* IORDRE : Maximum order of calculation of derivatives. */
/* NDIMEN : Dimension of the space. */
/* NCOEFF : Number of coefficients of the curve (degree+1). */
/* COURBE : Table of coefficients of the curve. */
/* TPARAM : Value of the parameter where the curve should be evaluated. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* TABPNT : Table of values of consecutive derivatives */
/* of parameter TPARAM. */
/* IERCOD : 0 = OK, */
/* 1 = incoherent input. */
/* COMMONS USED : */
/* ---------------- */
/* None. */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* The algorithm of calculation of derivatives is based on */
/* generalization of the Horner scheme : */
/* k 2 */
/* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
/* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
/* aj = a(j-1).x + u(k-j) */
/* bj = b(j-1).x + a(j-1) */
/* cj = c(j-1).x + b(j-1) */
/* So, it is obtained : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
/* The algorithm can be easily generalized for the calculation of */
/* (n) */
/* C (x) . */
/* --------- */
/* n! */
/* Reference : D. KNUTH, "The Art of Computer Programming" */
/* --------- Vol. 2/Seminumerical Algorithms */
/* Addison-Wesley Pub. Co. (1969) */
/* pages 423-425. */
/* ---> To evaluare derivatives at 0 and 1, it is preferable */
/* to use routine MDRV01.FOR . */
/* > */
/* **********************************************************************
*/
/* Name of the routine */
/* Parameter adjustments */
tabpnt_dim1 = *ndim;
--tabpnt;
courbe_dim1 = *ndim;
--courbe;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 2) {
AdvApp2Var_SysBase::mgenmsg_("MMDRVCB", 7L);
}
if (*ideriv < 0 || *ncoeff < 1) {
*iercod = 1;
goto L9999;
}
*iercod = 0;
/* ------------------- Initialization of table TABPNT -----------------
*/
ndgcrb = *ncoeff - 1;
i__1 = *ndim;
for (nd = 1; nd <= i__1; ++nd) {
tabpnt[nd] = courbe[nd + ndgcrb * courbe_dim1];
/* L100: */
}
if (*ideriv < 1) {
goto L200;
}
iptpnt = *ndim * *ideriv;
AdvApp2Var_SysBase::mvriraz_(&iptpnt,
&tabpnt[tabpnt_dim1 + 1]);
L200:
/* ------------------------ Calculation of parameter TPARAM ------------------
*/
i__1 = ndgcrb;
for (ndeg = 1; ndeg <= i__1; ++ndeg) {
i__2 = *ndim;
for (nd = 1; nd <= i__2; ++nd) {
for (i__ = *ideriv; i__ >= 1; --i__) {
tabpnt[nd + i__ * tabpnt_dim1] = tabpnt[nd + i__ *
tabpnt_dim1] * *tparam + tabpnt[nd + (i__ - 1) *
tabpnt_dim1];
/* L700: */
}
tabpnt[nd] = tabpnt[nd] * *tparam + courbe[nd + (ndgcrb - ndeg) *
courbe_dim1];
/* L600: */
}
/* L500: */
}
/* --------------------- Multiplication by factorial(I) -------------
*/
i__1 = *ideriv;
for (i__ = 2; i__ <= i__1; ++i__) {
i__2 = i__;
for (j = 2; j <= i__2; ++j) {
i__3 = *ndim;
for (nd = 1; nd <= i__3; ++nd) {
tabpnt[nd + i__ * tabpnt_dim1] = (doublereal) j * tabpnt[nd +
i__ * tabpnt_dim1];
/* L1200: */
}
/* L1100: */
}
/* L1000: */
}
/* --------------------------- The end ---------------------------------
*/
L9999:
if (*iercod > 0) {
AdvApp2Var_SysBase::maermsg_("MMDRVCB", iercod, 7L);
}
return 0;
} /* mmdrvcb_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmdrvck_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmdrvck_(integer *ncoeff,
integer *ndimen,
doublereal *courbe,
integer *ideriv,
doublereal *tparam,
doublereal *pntcrb)
{
/* Initialized data */
static doublereal mmfack[21] = { 1.,2.,6.,24.,120.,720.,5040.,40320.,
362880.,3628800.,39916800.,479001600.,6227020800.,87178291200.,
1.307674368e12,2.0922789888e13,3.55687428096e14,6.402373705728e15,
1.21645100408832e17,2.43290200817664e18,5.109094217170944e19 };
/* System generated locals */
integer courbe_dim1, courbe_offset, i__1, i__2;
/* Local variables */
integer i__, j, k, nd;
doublereal mfactk, bid;
/* IMPLICIT INTEGER (I-N) */
/* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
/* ***********************************************************************
*/
/* FONCTION : */
/* ---------- */
/* Calculate the value of a derived curve of order IDERIV in */
/* a point of parameter TPARAM. */
/* KEYWORDS : */
/* ----------- */
/* POSITIONING,CURVE,DERIVATIVE of ORDER K. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOEFF : Degree +1 of the curve. */
/* NDIMEN : Dimension of the space (2 or 3 in general) */
/* COURBE : Table of coefficients of the curve. */
/* IDERIV : Required order of derivation : 1=1st derivative, etc... */
/* TPARAM : Value of parameter of the curve. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* PNTCRB : Point of parameter TPARAM on the derivative of order */
/* IDERIV of CURVE. */
/* COMMONS USED : */
/* ---------------- */
/* MMCMCNP */
/* REFERENCES CALLED : */
/* ---------------------- */
/* None. */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* The code below was written basing on the following algorithm :
*/
/* Let P(t) = a1 + a2*t + ... an*t**n. The derivative of order k of P */
/* (containing n-k coefficients) is calculated as follows : */
/* Pk(t) = a(k+1)*CNP(k,k)*k! */
/* + a(k+2)*CNP(k+1,k)*k! * t */
/* . */
/* . */
/* . */
/* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
/* Evaluation is produced following the classic Horner scheme. */
/* > */
/* ***********************************************************************
*/
/* Factorials (1 to 21) caculated on VAX in R*16 */
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Serves to provide binomial coefficients (Pascal triangle). */
/* KEYWORDS : */
/* ----------- */
/* Binomial Coeff from 0 to 60. read only . init by block data */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* Binomial coefficients form a triangular matrix. */
/* This matrix is completed in table CNP by its transposition. */
/* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
/* Initialization is done by block-data MMLLL09.RES, */
/* created by program MQINICNP.FOR. */
/* > */
/* **********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
--pntcrb;
courbe_dim1 = *ndimen;
courbe_offset = courbe_dim1 + 1;
courbe -= courbe_offset;
/* Function Body */
/* -------------- Case when the order of derivative is greater than -------------------
*/
/* ---------------- the degree of the curve ---------------------
*/
if (*ideriv >= *ncoeff) {
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
pntcrb[nd] = 0.;
/* L100: */
}
goto L9999;
}
/* **********************************************************************
*/
/* General processing*/
/* **********************************************************************
*/
/* --------------------- Calculation of Factorial(IDERIV) ------------------
*/
k = *ideriv;
if (*ideriv <= 21 && *ideriv > 0) {
mfactk = mmfack[k - 1];
} else {
mfactk = 1.;
i__1 = k;
for (i__ = 2; i__ <= i__1; ++i__) {
mfactk *= i__;
/* L200: */
}
}
/* ------- Calculation of derivative of order IDERIV of CURVE in TPARAM -----
*/
/* ---> Attention : binomial coefficient C(n,m) is represented in */
/* MCCNP by CNP(N,M). */
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
pntcrb[nd] = courbe[nd + *ncoeff * courbe_dim1] * mmcmcnp_.cnp[*
ncoeff - 1 + k * 61] * mfactk;
/* L300: */
}
i__1 = k + 1;
for (j = *ncoeff - 1; j >= i__1; --j) {
bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
pntcrb[nd] = pntcrb[nd] * *tparam + courbe[nd + j * courbe_dim1] *
bid;
/* L500: */
}
/* L400: */
}
/* -------------------------------- The end -----------------------------
*/
L9999:
return 0 ;
} /* mmdrvck_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmeps1_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmeps1_(doublereal *epsilo)
{
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Extraction of EPS1 from COMMON MPRCSN. EPS1 is spatial zero */
/* equal to 1.D-9 */
/* KEYWORDS : */
/* ----------- */
/* MPRCSN,PRECISON,EPS1. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* None */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* EPSILO : Value of EPS1 (spatial zero (10**-9)) */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* EPS1 is ABSOLUTE spatial zero, so it is necessary */
/* to use it whenever it is necessary to test if a variable */
/* is null. For example, if the norm of a vector is lower than */
/* EPS1, this vector is NULL ! (when one works in */
/* REAL*8) It is absolutely not advised to test arguments */
/* compared to EPS1**2. Taking into account the rounding errors inevitable */
/* during calculations, this causes testing compared to 0.D0. */
/* > */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Gives tolerances of invalidity in stream */
/* as well as limits of iterative processes */
/* general context, modifiable by the user */
/* KEYWORDS : */
/* ----------- */
/* PARAMETER , TOLERANCE */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* INITIALISATION : profile , **VIA MPRFTX** at input in stream */
/* loading of default values of the profile in MPRFTX at input */
/* in stream. They are preserved in local variables of MPRFTX */
/* Reset of default values : MDFINT */
/* Interactive modification by the user : MDBINT */
/* ACCESS FUNCTION : MMEPS1 ... EPS1 */
/* MEPSPB ... EPS3,EPS4 */
/* MEPSLN ... EPS2, NITERM , NITERR */
/* MEPSNR ... EPS2 , NITERM */
/* MITERR ... NITERR */
/* > */
/* ***********************************************************************
*/
/* NITERM : max nb of iterations */
/* NITERR : nb of rapid iterations */
/* EPS1 : tolerance of 3D null distance */
/* EPS2 : tolerance of parametric null distance */
/* EPS3 : tolerance to avoid division by 0.. */
/* EPS4 : angular tolerance */
/* ***********************************************************************
*/
*epsilo = mmprcsn_.eps1;
return 0 ;
} /* mmeps1_ */
//=======================================================================
//function : mmexthi_
//purpose :
//=======================================================================
int mmexthi_(integer *ndegre,
doublereal *hwgaus)
{
/* System generated locals */
integer i__1;
/* Local variables */
integer iadd, ideb, ndeg2, nmod2, ii, ibb;
integer kpt;
/* **********************************************************************
*/
/* FONCTION : */
/* ---------- */
/* Extract of common LDGRTL the weight of formulas of */
/* Gauss quadrature on all roots of Legendre polynoms of degree */
/* NDEGRE defined on [-1,1]. */
/* KEYWORDS : */
/* ----------- */
/* ALL, AB_SPECIFI::COMMON&, EXTRACTION, &WEIGHT, &GAUSS. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDEGRE : Mathematic degree of Legendre polynom. It should have */
/* 2 <= NDEGRE <= 61. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* HWGAUS : The table of weights of Gauss quadrature formulas */
/* relative to NDEGRE roots of a polynome de Legendre de */
/* degre NDEGRE. */
/* COMMONS UTILISES : */
/* ---------------- */
/* MLGDRTL */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ATTENTION: The condition on NDEGRE ( 2 <= NDEGRE <= 61) is not */
/* tested. The caller should make the test. */
/* Name of the routine */
/* Common MLGDRTL: */
/* This common includes POSITIVE roots of Legendre polynims */
/* AND weights of Gauss quadrature formulas on all */
/* POSITIVE roots of Legendre polynoms. */
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* The common of Legendre roots. */
/* KEYWORDS : */
/* ----------- */
/* BASE LEGENDRE */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* ROOTAB : Table of all roots of Legendre polynoms */
/* within the interval [0,1]. They are ranked for the degrees increasing from */
/* 2 to 61. */
/* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
/* The adressing is the same. */
/* HI0TAB : Table of Legendre interpolators for root x=0 */
/* of polynoms of UNEVEN degree. */
/* RTLTB0 : Table of Li(uk) where uk are the roots of */
/* Legendre polynom of EVEN degree. */
/* RTLTB1 : Table of Li(uk) where uk are the roots of */
/* Legendre polynom of UNEVEN degree. */
/************************************************************************
*****/
/* Parameter adjustments */
--hwgaus;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 3) {
AdvApp2Var_SysBase::mgenmsg_("MMEXTHI", 7L);
}
ndeg2 = *ndegre / 2;
nmod2 = *ndegre % 2;
/* Address of Gauss weight associated to the 1st strictly */
/* positive root of Legendre polynom of degree NDEGRE in MLGDRTL. */
iadd = ndeg2 * (ndeg2 - 1) / 2 + 1;
/* Index of the 1st HWGAUS element associated to the 1st strictly */
/* positive root of Legendre polynom of degree NDEGRE. */
ideb = (*ndegre + 1) / 2 + 1;
/* Reading of weights associated to strictly positive roots. */
i__1 = *ndegre;
for (ii = ideb; ii <= i__1; ++ii) {
kpt = iadd + ii - ideb;
hwgaus[ii] = mlgdrtl_.hiltab[kpt + nmod2 * 465 - 1];
/* L100: */
}
/* For strictly negative roots, the weight is the same. */
/* i.e HW(1) = HW(NDEGRE), HW(2) = HW(NDEGRE-1), etc... */
i__1 = ndeg2;
for (ii = 1; ii <= i__1; ++ii) {
hwgaus[ii] = hwgaus[*ndegre + 1 - ii];
/* L200: */
}
/* Case of uneven NDEGRE, 0 is root of Legendre polynom, */
/* associated Gauss weights are loaded. */
if (nmod2 == 1) {
hwgaus[ndeg2 + 1] = mlgdrtl_.hi0tab[ndeg2];
}
/* --------------------------- The end ----------------------------------
*/
if (ibb >= 3) {
AdvApp2Var_SysBase::mgsomsg_("MMEXTHI", 7L);
}
return 0;
} /* mmexthi_ */
//=======================================================================
//function : mmextrl_
//purpose :
//=======================================================================
int mmextrl_(integer *ndegre,
doublereal *rootlg)
{
/* System generated locals */
integer i__1;
/* Local variables */
integer iadd, ideb, ndeg2, nmod2, ii, ibb;
integer kpt;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Extract of the Common LDGRTL of Legendre polynom roots */
/* of degree NDEGRE defined on [-1,1]. */
/* KEYWORDS : */
/* ----------- */
/* ALL, AB_SPECIFI::COMMON&, EXTRACTION, &ROOT, &LEGENDRE. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDEGRE : Mathematic degree of Legendre polynom. */
/* It is required to have 2 <= NDEGRE <= 61. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* ROOTLG : The table of roots of Legendre polynom of degree */
/* NDEGRE defined on [-1,1]. */
/* COMMONS USED : */
/* ---------------- */
/* MLGDRTL */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ATTENTION: Condition of NDEGRE ( 2 <= NDEGRE <= 61) is not */
/* tested. The caller should make the test. */
/* > */
/* **********************************************************************
*/
/* Name of the routine */
/* Common MLGDRTL: */
/* This common includes POSITIVE roots of Legendre polynoms */
/* AND the weight of Gauss quadrature formulas on all */
/* POSITIVE roots of Legendre polynoms. */
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* The common of Legendre roots. */
/* KEYWORDS : */
/* ----------- */
/* BASE LEGENDRE */
/* ***********************************************************************
*/
/* ROOTAB : Table of all roots of Legendre polynoms */
/* within the interval [0,1]. They are ranked for the degrees increasing from */
/* 2 to 61. */
/* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
/* The adressing is the same. */
/* HI0TAB : Table of Legendre interpolators for root x=0 */
/* of polynoms of UNEVEN degree. */
/* RTLTB0 : Table of Li(uk) where uk are the roots of */
/* Legendre polynom of EVEN degree. */
/* RTLTB1 : Table of Li(uk) where uk are the roots of */
/* Legendre polynom of UNEVEN degree. */
/************************************************************************
*****/
/* Parameter adjustments */
--rootlg;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 3) {
AdvApp2Var_SysBase::mgenmsg_("MMEXTRL", 7L);
}
ndeg2 = *ndegre / 2;
nmod2 = *ndegre % 2;
/* Address of the 1st strictly positive root of Legendre polynom */
/* of degree NDEGRE in MLGDRTL. */
iadd = ndeg2 * (ndeg2 - 1) / 2 + 1;
/* Indice, in ROOTLG, of the 1st strictly positive root */
/* of Legendre polynom of degree NDEGRE. */
ideb = (*ndegre + 1) / 2 + 1;
/* Reading of strictly positive roots. */
i__1 = *ndegre;
for (ii = ideb; ii <= i__1; ++ii) {
kpt = iadd + ii - ideb;
rootlg[ii] = mlgdrtl_.rootab[kpt + nmod2 * 465 - 1];
/* L100: */
}
/* Strictly negative roots are equal to positive roots
*/
/* to the sign i.e RT(1) = -RT(NDEGRE), RT(2) = -RT(NDEGRE-1), etc...
*/
i__1 = ndeg2;
for (ii = 1; ii <= i__1; ++ii) {
rootlg[ii] = -rootlg[*ndegre + 1 - ii];
/* L200: */
}
/* Case NDEGRE uneven, 0 is root of Legendre polynom. */
if (nmod2 == 1) {
rootlg[ndeg2 + 1] = 0.;
}
/* -------------------------------- THE END -----------------------------
*/
if (ibb >= 3) {
AdvApp2Var_SysBase::mgenmsg_("MMEXTRL", 7L);
}
return 0;
} /* mmextrl_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmfmca8_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmfmca8_(const integer *ndimen,
const integer *ncoefu,
const integer *ncoefv,
const integer *ndimax,
const integer *ncfumx,
const integer *,//ncfvmx,
doublereal *tabini,
doublereal *tabres)
{
/* System generated locals */
integer tabini_dim1, tabini_dim2, tabini_offset, tabres_dim1, tabres_dim2,
tabres_offset;
/* Local variables */
integer i__, j, k, ilong;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Expansion of a table containing only most important things into a */
/* greater data table. */
/* KEYWORDS : */
/* ----------- */
/* ALL, MATH_ACCES:: CARREAU&, DECOMPRESSION, &CARREAU */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMEN: Dimension of the workspace. */
/* NCOEFU: Degree +1 of the table by u. */
/* NCOEFV: Degree +1 of the table by v. */
/* NDIMAX: Max dimension of the space. */
/* NCFUMX: Max Degree +1 of the table by u. */
/* NCFVMX: Max Degree +1 of the table by v. */
/* TABINI: The table to be decompressed. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* TABRES: Decompressed table. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* The following call : */
/* CALL MMFMCA8(NDIMEN,NCOEFU,NCOEFV,NDIMAX,NCFUMX,NCFVMX,TABINI,TABINI)
*/
/* where TABINI is input/output argument, is possible provided */
/* that the caller has declared TABINI in (NDIMAX,NCFUMX,NCFVMX) */
/* ATTENTION : it is not checked that NDIMAX >= NDIMEN, */
/* NCOEFU >= NCFMXU and NCOEFV >= NCFMXV. */
/* > */
/* **********************************************************************
*/
/* Parameter adjustments */
tabini_dim1 = *ndimen;
tabini_dim2 = *ncoefu;
tabini_offset = tabini_dim1 * (tabini_dim2 + 1) + 1;
tabini -= tabini_offset;
tabres_dim1 = *ndimax;
tabres_dim2 = *ncfumx;
tabres_offset = tabres_dim1 * (tabres_dim2 + 1) + 1;
tabres -= tabres_offset;
/* Function Body */
if (*ndimax == *ndimen) {
goto L1000;
}
/* ----------------------- decompression NDIMAX<>NDIMEN -----------------
*/
for (k = *ncoefv; k >= 1; --k) {
for (j = *ncoefu; j >= 1; --j) {
for (i__ = *ndimen; i__ >= 1; --i__) {
tabres[i__ + (j + k * tabres_dim2) * tabres_dim1] = tabini[
i__ + (j + k * tabini_dim2) * tabini_dim1];
/* L300: */
}
/* L200: */
}
/* L100: */
}
goto L9999;
/* ----------------------- decompression NDIMAX=NDIMEN ------------------
*/
L1000:
if (*ncoefu == *ncfumx) {
goto L2000;
}
ilong = (*ndimen << 3) * *ncoefu;
for (k = *ncoefv; k >= 1; --k) {
AdvApp2Var_SysBase::mcrfill_(&ilong,
&tabini[(k * tabini_dim2 + 1) * tabini_dim1 + 1],
&tabres[(k * tabres_dim2 + 1) * tabres_dim1 + 1]);
/* L500: */
}
goto L9999;
/* ----------------- decompression NDIMAX=NDIMEN,NCOEFU=NCFUMX ----------
*/
L2000:
ilong = (*ndimen << 3) * *ncoefu * *ncoefv;
AdvApp2Var_SysBase::mcrfill_(&ilong,
&tabini[tabini_offset],
&tabres[tabres_offset]);
goto L9999;
/* ---------------------------- The end ---------------------------------
*/
L9999:
return 0;
} /* mmfmca8_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmfmca9_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmfmca9_(integer *ndimax,
integer *ncfumx,
integer *,//ncfvmx,
integer *ndimen,
integer *ncoefu,
integer *ncoefv,
doublereal *tabini,
doublereal *tabres)
{
/* System generated locals */
integer tabini_dim1, tabini_dim2, tabini_offset, tabres_dim1, tabres_dim2,
tabres_offset, i__1, i__2, i__3;
/* Local variables */
integer i__, j, k, ilong;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Compression of a data table in a table */
/* containing only the main data (the input table is not removed). */
/* KEYWORDS: */
/* ----------- */
/* ALL, MATH_ACCES:: CARREAU&, COMPRESSION, &CARREAU */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMAX: Max dimension of the space. */
/* NCFUMX: Max degree +1 of the table by u. */
/* NCFVMX: Max degree +1 of the table by v. */
/* NDIMEN: Dimension of the workspace. */
/* NCOEFU: Degree +1 of the table by u. */
/* NCOEFV: Degree +1 of the table by v. */
/* TABINI: The table to compress. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* TABRES: The compressed table. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* The following call : */
/* CALL MMFMCA9(NDIMAX,NCFUMX,NCFVMX,NDIMEN,NCOEFU,NCOEFV,TABINI,TABINI)
*/
/* where TABINI is input/output argument, is possible provided */
/* that the caller has checked that : */
/* NDIMAX > NDIMEN, */
/* or NDIMAX = NDIMEN and NCFUMX > NCOEFU */
/* or NDIMAX = NDIMEN, NCFUMX = NCOEFU and NCFVMX > NCOEFV */
/* These conditions are not tested in the program. */
/* > */
/* **********************************************************************
*/
/* Parameter adjustments */
tabini_dim1 = *ndimax;
tabini_dim2 = *ncfumx;
tabini_offset = tabini_dim1 * (tabini_dim2 + 1) + 1;
tabini -= tabini_offset;
tabres_dim1 = *ndimen;
tabres_dim2 = *ncoefu;
tabres_offset = tabres_dim1 * (tabres_dim2 + 1) + 1;
tabres -= tabres_offset;
/* Function Body */
if (*ndimen == *ndimax) {
goto L1000;
}
/* ----------------------- Compression NDIMEN<>NDIMAX -------------------
*/
i__1 = *ncoefv;
for (k = 1; k <= i__1; ++k) {
i__2 = *ncoefu;
for (j = 1; j <= i__2; ++j) {
i__3 = *ndimen;
for (i__ = 1; i__ <= i__3; ++i__) {
tabres[i__ + (j + k * tabres_dim2) * tabres_dim1] = tabini[
i__ + (j + k * tabini_dim2) * tabini_dim1];
/* L300: */
}
/* L200: */
}
/* L100: */
}
goto L9999;
/* ----------------------- Compression NDIMEN=NDIMAX --------------------
*/
L1000:
if (*ncoefu == *ncfumx) {
goto L2000;
}
ilong = (*ndimen << 3) * *ncoefu;
i__1 = *ncoefv;
for (k = 1; k <= i__1; ++k) {
AdvApp2Var_SysBase::mcrfill_(&ilong,
&tabini[(k * tabini_dim2 + 1) * tabini_dim1 + 1],
&tabres[(k * tabres_dim2 + 1) * tabres_dim1 + 1]);
/* L500: */
}
goto L9999;
/* ----------------- Compression NDIMEN=NDIMAX,NCOEFU=NCFUMX ------------
*/
L2000:
ilong = (*ndimen << 3) * *ncoefu * *ncoefv;
AdvApp2Var_SysBase::mcrfill_(&ilong,
&tabini[tabini_offset],
&tabres[tabres_offset]);
goto L9999;
/* ---------------------------- The end ---------------------------------
*/
L9999:
return 0;
} /* mmfmca9_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmfmcar_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmfmcar_(integer *ndimen,
integer *ncofmx,
integer *ncoefu,
integer *ncoefv,
doublereal *patold,
doublereal *upara1,
doublereal *upara2,
doublereal *vpara1,
doublereal *vpara2,
doublereal *patnew,
integer *iercod)
{
integer c__8 = 8;
/* System generated locals */
integer patold_dim1, patold_dim2, patnew_dim1, patnew_dim2,
i__1, patold_offset,patnew_offset;
/* Local variables */
doublereal* tbaux = 0;
integer ksize, numax, kk;
intptr_t iofst;
integer ibb, ier;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* LIMITATION OF A SQUARE DEFINED ON (0,1)*(0,1) BETWEEN ISOS */
/* UPARA1 AND UPARA2 (BY U) AND VPARA1 AND VPARA2 BY V. */
/* KEYWORDS : */
/* ----------- */
/* LIMITATION , SQUARE , PARAMETER */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX: MAX NUMBER OF COEFF OF THE SQUARE BY U */
/* NCOEFU: NUMBER OF COEFF OF THE SQUARE BY U */
/* NCOEFV: NUMBER OF COEFF OF THE SQUARE BY V */
/* PATOLD : THE SQUARE IS LIMITED BY UPARA1,UPARA2 AND VPARA1,VPARA2
.*/
/* UPARA1 : LOWER LIMIT OF U */
/* UPARA2 : UPPER LIMIT OF U */
/* VPARA1 : LOWER LIMIT OF V */
/* VPARA2 : UPPER LIMIT OF V */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* PATNEW : RELIMITED SQUARE, DEFINED ON (0,1)**2 */
/* IERCOD : =10 COEFF NB TOO GREAT OR NULL */
/* =13 PB IN THE DYNAMIC ALLOCATION */
/* = 0 OK. */
/* COMMONS USED : */
/* ---------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ---> The following call : */
/* CALL MMFMCAR(NCOFMX,NCOEFU,NCOEFV,PATOLD,UPARA1,UPARA2,VPARA1,VPARA2
*/
/* ,PATOLD), */
/* where PATOLD is input/output argument is absolutely legal. */
/* ---> The max number of coeff by u and v of PATOLD is 61 */
/* ---> If NCOEFU < NCOFMX, the data is compressed by MMFMCA9 before */
/* limitation by v to get time during the execution */
/* of MMARC41 that follows (the square is processed as a curve of
*/
/* dimension NDIMEN*NCOEFU possessing NCOEFV coefficients). */
/* > */
/* ***********************************************************************
*/
/* Name of the routine */
/* Parameter adjustments */
patnew_dim1 = *ndimen;
patnew_dim2 = *ncofmx;
patnew_offset = patnew_dim1 * (patnew_dim2 + 1) + 1;
patnew -= patnew_offset;
patold_dim1 = *ndimen;
patold_dim2 = *ncofmx;
patold_offset = patold_dim1 * (patold_dim2 + 1) + 1;
patold -= patold_offset;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 2) {
AdvApp2Var_SysBase::mgenmsg_("MMFMCAR", 7L);
}
*iercod = 0;
iofst = 0;
AdvApp2Var_SysBase anAdvApp2Var_SysBase;
/* **********************************************************************
*/
/* TEST OF COEFFICIENT NUMBERS */
/* **********************************************************************
*/
if (*ncofmx < *ncoefu) {
*iercod = 10;
goto L9999;
}
if (*ncoefu < 1 || *ncoefu > 61 || *ncoefv < 1 || *ncoefv > 61) {
*iercod = 10;
goto L9999;
}
/* **********************************************************************
*/
/* CASE WHEN UPARA1=VPARA1=0 AND UPARA2=VPARA2=1 */
/* **********************************************************************
*/
if (*upara1 == 0. && *upara2 == 1. && *vpara1 == 0. && *vpara2 == 1.) {
ksize = (*ndimen << 3) * *ncofmx * *ncoefv;
AdvApp2Var_SysBase::mcrfill_(&ksize,
&patold[patold_offset],
&patnew[patnew_offset]);
goto L9999;
}
/* **********************************************************************
*/
/* LIMITATION BY U */
/* **********************************************************************
*/
if (*upara1 == 0. && *upara2 == 1.) {
goto L2000;
}
i__1 = *ncoefv;
for (kk = 1; kk <= i__1; ++kk) {
mmarc41_(ndimen, ndimen, ncoefu, &patold[(kk * patold_dim2 + 1) *
patold_dim1 + 1], upara1, upara2, &patnew[(kk * patnew_dim2 +
1) * patnew_dim1 + 1], iercod);
/* L100: */
}
/* **********************************************************************
*/
/* LIMITATION BY V */
/* **********************************************************************
*/
L2000:
if (*vpara1 == 0. && *vpara2 == 1.) {
goto L9999;
}
/* ----------- LIMITATION BY V (WITH COMPRESSION I.E. NCOEFU<NCOFMX) ----
*/
numax = *ndimen * *ncoefu;
if (*ncofmx != *ncoefu) {
/* ------------------------- Dynamic allocation -------------------
---- */
ksize = *ndimen * *ncoefu * *ncoefv;
anAdvApp2Var_SysBase.mcrrqst_(&c__8, &ksize, tbaux, &iofst, &ier);
if (ier > 0) {
*iercod = 13;
goto L9900;
}
/* --------------- Compression by (NDIMEN,NCOEFU,NCOEFV) ------------
---- */
if (*upara1 == 0. && *upara2 == 1.) {
AdvApp2Var_MathBase::mmfmca9_(ndimen,
ncofmx,
ncoefv,
ndimen,
ncoefu,
ncoefv,
&patold[patold_offset],
&tbaux[iofst]);
} else {
AdvApp2Var_MathBase::mmfmca9_(ndimen,
ncofmx,
ncoefv,
ndimen,
ncoefu,
ncoefv,
&patnew[patnew_offset],
&tbaux[iofst]);
}
/* ------------------------- Limitation by v ------------------------
---- */
mmarc41_(&numax, &numax, ncoefv, &tbaux[iofst], vpara1, vpara2, &
tbaux[iofst], iercod);
/* --------------------- Expansion of TBAUX into PATNEW -------------
--- */
AdvApp2Var_MathBase::mmfmca8_(ndimen, ncoefu, ncoefv, ndimen, ncofmx, ncoefv, &tbaux[iofst]
, &patnew[patnew_offset]);
goto L9900;
/* -------- LIMITATION BY V (WITHOUT COMPRESSION I.E. NCOEFU=NCOFMX) ---
---- */
} else {
if (*upara1 == 0. && *upara2 == 1.) {
mmarc41_(&numax, &numax, ncoefv, &patold[patold_offset], vpara1,
vpara2, &patnew[patnew_offset], iercod);
} else {
mmarc41_(&numax, &numax, ncoefv, &patnew[patnew_offset], vpara1,
vpara2, &patnew[patnew_offset], iercod);
}
goto L9999;
}
/* **********************************************************************
*/
/* DESALLOCATION */
/* **********************************************************************
*/
L9900:
if (iofst != 0) {
anAdvApp2Var_SysBase.mcrdelt_(&c__8, &ksize, tbaux, &iofst, &ier);
}
if (ier > 0) {
*iercod = 13;
}
/* ------------------------------ The end -------------------------------
*/
L9999:
if (*iercod > 0) {
AdvApp2Var_SysBase::maermsg_("MMFMCAR", iercod, 7L);
}
if (ibb >= 2) {
AdvApp2Var_SysBase::mgsomsg_("MMFMCAR", 7L);
}
return 0;
} /* mmfmcar_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmfmcb5_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmfmcb5_(integer *isenmsc,
integer *ndimax,
integer *ncf1mx,
doublereal *courb1,
integer *ncoeff,
integer *ncf2mx,
integer *ndimen,
doublereal *courb2,
integer *iercod)
{
/* System generated locals */
integer courb1_dim1, courb1_offset, courb2_dim1, courb2_offset, i__1,
i__2;
/* Local variables */
integer i__, nboct, nd;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Reformating (and eventual compression/decompression) of curve */
/* (ndim,.) by (.,ndim) and vice versa. */
/* KEYWORDS : */
/* ----------- */
/* ALL , MATH_ACCES :: */
/* COURBE&, REORGANISATION,COMPRESSION,INVERSION , &COURBE */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* ISENMSC : required direction of the transfer : */
/* 1 : passage of (NDIMEN,.) ---> (.,NDIMEN) direction to AB
*/
/* -1 : passage of (.,NDIMEN) ---> (NDIMEN,.) direction to TS,T
V*/
/* NDIMAX : format / dimension */
/* NCF1MX : format by t of COURB1 */
/* if ISENMSC= 1 : COURB1: The curve to be processed (NDIMAX,.) */
/* NCOEFF : number of coeff of the curve */
/* NCF2MX : format by t of COURB2 */
/* NDIMEN : dimension of the curve and format of COURB2 */
/* if ISENMSC=-1 : COURB2: The curve to be processed (.,NDIMEN) */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* if ISENMSC= 1 : COURB2: The resulting curve (.,NDIMEN) */
/* if ISENMSC=-1 : COURB1: The resulting curve (NDIMAX,.) */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* allow to process the usual transfers as follows : */
/* | ---- ISENMSC = 1 ---- | | ---- ISENMSC =-1 ----- | */
/* TS (3,21) --> (21,3) AB ; AB (21,3) --> (3,21) TS */
/* TS (3,21) --> (NU,3) AB ; AB (NU,3) --> (3,21) TS */
/* (3,NU) --> (21,3) AB ; AB (21,3) --> (3,NU) */
/* (3,NU) --> (NU,3) AB ; AB (NU,3) --> (3,NU) */
/* > */
/* ***********************************************************************
*/
/* Parameter adjustments */
courb1_dim1 = *ndimax;
courb1_offset = courb1_dim1 + 1;
courb1 -= courb1_offset;
courb2_dim1 = *ncf2mx;
courb2_offset = courb2_dim1 + 1;
courb2 -= courb2_offset;
/* Function Body */
if (*ndimen > *ndimax || *ncoeff > *ncf1mx || *ncoeff > *ncf2mx) {
goto L9119;
}
if (*ndimen == 1 && *ncf1mx == *ncf2mx) {
nboct = *ncf2mx << 3;
if (*isenmsc == 1) {
AdvApp2Var_SysBase::mcrfill_(&nboct,
&courb1[courb1_offset],
&courb2[courb2_offset]);
}
if (*isenmsc == -1) {
AdvApp2Var_SysBase::mcrfill_(&nboct,
&courb2[courb2_offset],
&courb1[courb1_offset]);
}
*iercod = -3136;
goto L9999;
}
*iercod = 0;
if (*isenmsc == 1) {
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
i__2 = *ncoeff;
for (i__ = 1; i__ <= i__2; ++i__) {
courb2[i__ + nd * courb2_dim1] = courb1[nd + i__ *
courb1_dim1];
/* L400: */
}
/* L500: */
}
} else if (*isenmsc == -1) {
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
i__2 = *ncoeff;
for (i__ = 1; i__ <= i__2; ++i__) {
courb1[nd + i__ * courb1_dim1] = courb2[i__ + nd *
courb2_dim1];
/* L1400: */
}
/* L1500: */
}
} else {
*iercod = 3164;
}
goto L9999;
/* ***********************************************************************
*/
L9119:
*iercod = 3119;
L9999:
if (*iercod != 0) {
AdvApp2Var_SysBase::maermsg_("MMFMCB5", iercod, 7L);
}
return 0;
} /* mmfmcb5_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmfmtb1_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmfmtb1_(integer *maxsz1,
doublereal *table1,
integer *isize1,
integer *jsize1,
integer *maxsz2,
doublereal *table2,
integer *isize2,
integer *jsize2,
integer *iercod)
{
integer c__8 = 8;
/* System generated locals */
integer table1_dim1, table1_offset, table2_dim1, table2_offset, i__1,
i__2;
/* Local variables */
doublereal* work = 0;
integer ilong, isize, ii, jj, ier = 0;
intptr_t iofst = 0,iipt, jjpt;
/************************************************************************
*******/
/* FUNCTION : */
/* ---------- */
/* Inversion of elements of a rectangular table (T1(i,j) */
/* loaded in T2(j,i)) */
/* KEYWORDS : */
/* ----------- */
/* ALL, MATH_ACCES :: TABLEAU&, INVERSION, &TABLEAU */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* MAXSZ1: Max Nb of elements by the 1st dimension of TABLE1. */
/* TABLE1: Table of reals by two dimensions. */
/* ISIZE1: Nb of useful elements of TABLE1 on the 1st dimension */
/* JSIZE1: Nb of useful elements of TABLE1 on the 2nd dimension */
/* MAXSZ2: Nb max of elements by the 1st dimension of TABLE2. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* TABLE2: Table of reals by two dimensions, containing the transposition */
/* of the rectangular table TABLE1. */
/* ISIZE2: Nb of useful elements of TABLE2 on the 1st dimension */
/* JSIZE2: Nb of useful elements of TABLE2 on the 2nd dimension */
/* IERCOD: Erroe coder. */
/* = 0, ok. */
/* = 1, error in the dimension of tables */
/* ether MAXSZ1 < ISIZE1 (table TABLE1 too small). */
/* or MAXSZ2 < JSIZE1 (table TABLE2 too small). */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ---------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* It is possible to use TABLE1 as input and output table i.e. */
/* call: */
/* CALL MMFMTB1(MAXSZ1,TABLE1,ISIZE1,JSIZE1,MAXSZ2,TABLE1 */
/* ,ISIZE2,JSIZE2,IERCOD) */
/* is valuable. */
/* > */
/* **********************************************************************
*/
/* Parameter adjustments */
table1_dim1 = *maxsz1;
table1_offset = table1_dim1 + 1;
table1 -= table1_offset;
table2_dim1 = *maxsz2;
table2_offset = table2_dim1 + 1;
table2 -= table2_offset;
AdvApp2Var_SysBase anAdvApp2Var_SysBase;
/* Function Body */
*iercod = 0;
if (*isize1 > *maxsz1 || *jsize1 > *maxsz2) {
goto L9100;
}
iofst = 0;
isize = *maxsz2 * *isize1;
anAdvApp2Var_SysBase.mcrrqst_(&c__8, &isize, work, &iofst, &ier);
if (ier > 0) {
goto L9200;
}
/* DO NOT BE AFRAID OF CRUSHING. */
i__1 = *isize1;
for (ii = 1; ii <= i__1; ++ii) {
iipt = (ii - 1) * *maxsz2 + iofst;
i__2 = *jsize1;
for (jj = 1; jj <= i__2; ++jj) {
jjpt = iipt + (jj - 1);
work[jjpt] = table1[ii + jj * table1_dim1];
/* L200: */
}
/* L100: */
}
ilong = isize << 3;
AdvApp2Var_SysBase::mcrfill_(&ilong,
&work[iofst],
&table2[table2_offset]);
/* -------------- The number of elements of TABLE2 is returned ------------
*/
ii = *isize1;
*isize2 = *jsize1;
*jsize2 = ii;
goto L9999;
/* ------------------------------- THE END ------------------------------
*/
/* --> Invalid input. */
L9100:
*iercod = 1;
goto L9999;
/* --> Pb of allocation. */
L9200:
*iercod = 2;
goto L9999;
L9999:
if (iofst != 0) {
anAdvApp2Var_SysBase.mcrdelt_(&c__8, &isize, work, &iofst, &ier);
}
if (ier > 0) {
*iercod = 2;
}
return 0;
} /* mmfmtb1_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmgaus1_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmgaus1_(integer *ndimf,
int (*bfunx) (
integer *ninteg,
doublereal *parame,
doublereal *vfunj1,
integer *iercod
),
integer *k,
doublereal *xd,
doublereal *xf,
doublereal *saux1,
doublereal *saux2,
doublereal *somme,
integer *niter,
integer *iercod)
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
integer ndeg;
doublereal h__[20];
integer j;
doublereal t, u[20], x;
integer idimf;
doublereal c1x, c2x;
/* **********************************************************************
*/
/* FUNCTION : */
/* -------- */
/* Calculate the integral of function BFUNX passed in parameter */
/* between limits XD and XF . */
/* The function should be calculated for any value */
/* of the variable in the given interval.. */
/* The method GAUSS-LEGENDRE is used. */
/* For explications refer to the book : */
/* Complements de mathematiques a l'usage des Ingenieurs de */
/* l'electrotechnique et des telecommunications. */
/* Par Andre ANGOT - Collection technique et scientifique du CNET
*/
/* page 772 .... */
/* The degree of LEGENDRE polynoms used is passed in parameter.
*/
/* KEYWORDS : */
/* --------- */
/* INTEGRATION,LEGENDRE,GAUSS */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMF : Dimension of the function */
/* BFUNX : Function to integrate passed as argument */
/* Should be declared as EXTERNAL in the call routine. */
/* SUBROUTINE BFUNX(NDIMF,X,VAL,IER) */
/* REAL *8 X,VAL */
/* K : Parameter determining the degree of the LEGENDRE polynom that
*/
/* can take a value between 0 and 10. */
/* The degree of the polynom is equal to 4 k, that is 4, 8,
*/
/* 12, 16, 20, 24, 28, 32, 36 and 40. */
/* If K is not correct, the degree is set to 40 directly.
*/
/* XD : Lower limit of the interval of integration. */
/* XF : Upper limit of the interval of integration. */
/* SAUX1 : Auxiliary table */
/* SAUX2 : Auxiliary table */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* SOMME : Value of the integral */
/* NITER : Number of iterations to be carried out. */
/* It is equal to the degree of the polynom. */
/* IER : Error code : */
/* < 0 ==> Attention - Warning */
/* = 0 ==> Everything is OK */
/* > 0 ==> Critical error - Apply special processing */
/* ==> Error in the calculation of BFUNX (return code */
/* of this routine */
/* If error => SUM = 0 */
/* COMMONS USED : */
/* ----------------- */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Type Name */
/* @ BFUNX MVGAUS0 */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* --------------------------------- */
/* See the explanations detailed in the listing */
/* Use of the GAUSS method (orthogonal polynoms) */
/* The symmetry of roots of these polynomes is used */
/* Depending on K, the degree of the interpolated polynom grows.
*/
/* If you wish to calculate the integral with a given precision, */
/* loop on k varying from 1 to 10 and test the difference of 2
*/
/* consecutive iterations. Stop the loop if this difference is less that */
/* an epsilon value set to 10E-6 for example. */
/* If S1 and S2 are 2 successive iterations, test following this example :
*/
/* AF=DABS(S1-S2) */
/* AS=DABS(S2) */
/* If AS < 1 test if FS < eps otherwise test if AF/AS < eps
*/
/* -- ----- ----- */
/* > */
/************************************************************************
******/
/* DECLARATIONS */
/************************************************************************
******/
/* ****** General Initialization */
/* Parameter adjustments */
--somme;
--saux2;
--saux1;
/* Function Body */
AdvApp2Var_SysBase::mvriraz_(ndimf,
&somme[1]);
*iercod = 0;
/* ****** Loading of coefficients U and H ** */
/* -------------------------------------------- */
mvgaus0_(k, u, h__, &ndeg, iercod);
if (*iercod > 0) {
goto L9999;
}
/* ****** C1X => Medium interval point [XD,XF] */
/* ****** C2X => 1/2 amplitude interval [XD,XF] */
c1x = (*xf + *xd) * .5;
c2x = (*xf - *xd) * .5;
/* ---------------------------------------- */
/* ****** Integration for degree NDEG ** */
/* ---------------------------------------- */
i__1 = ndeg;
for (j = 1; j <= i__1; ++j) {
t = c2x * u[j - 1];
x = c1x + t;
(*bfunx)(ndimf, &x, &saux1[1], iercod);
if (*iercod != 0) {
goto L9999;
}
x = c1x - t;
(*bfunx)(ndimf, &x, &saux2[1], iercod);
if (*iercod != 0) {
goto L9999;
}
i__2 = *ndimf;
for (idimf = 1; idimf <= i__2; ++idimf) {
somme[idimf] += h__[j - 1] * (saux1[idimf] + saux2[idimf]);
}
}
*niter = ndeg << 1;
i__1 = *ndimf;
for (idimf = 1; idimf <= i__1; ++idimf) {
somme[idimf] *= c2x;
}
/* ****** End of sub-program ** */
L9999:
return 0 ;
} /* mmgaus1_ */
//=======================================================================
//function : mmherm0_
//purpose :
//=======================================================================
int mmherm0_(doublereal *debfin,
integer *iercod)
{
integer c__576 = 576;
integer c__6 = 6;
/* System generated locals */
integer i__1, i__2;
doublereal d__1;
/* Local variables */
doublereal amat[36] /* was [6][6] */;
integer iord[2];
doublereal prod;
integer iord1, iord2;
doublereal miden[36] /* was [6][6] */;
integer ncmat;
doublereal epspi, d1, d2;
integer ii, jj, pp, ncf;
doublereal cof[6];
integer iof[2], ier;
doublereal mat[36] /* was [6][6] */;
integer cot;
doublereal abid[72] /* was [12][6] */;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* INIT OF COEFFS. OF POLYNOMS OF HERMIT INTERPOLATION */
/* KEYWORDS : */
/* ----------- */
/* MATH_ACCES :: HERMITE */
/* INPUT ARGUMENTS */
/* -------------------- */
/* DEBFIN : PARAMETERS DEFINING THE CONSTRAINTS */
/* DEBFIN(1) : FIRST PARAMETER */
/* DEBFIN(2) : SECOND PARAMETER */
/* ONE SHOULD HAVE: */
/* ABS (DEBFIN(I)) < 100 */
/* and */
/* (ABS(DEBFIN(1)+ABS(DEBFIN(2))) > 1/100 */
/* (for overflows) */
/* ABS(DEBFIN(2)-DEBFIN(1)) / (ABS(DEBFIN(1)+ABS(DEBFIN(2))) > 1/100
*/
/* (for the conditioning) */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* IERCOD : Error code : 0 : O.K. */
/* 1 : value of DEBFIN */
/* are unreasonable */
/* -1 : init was already done */
/* (OK but no processing) */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Type Name */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* This program initializes the coefficients of Hermit polynoms */
/* that are read later by MMHERM1 */
/* ***********************************************************************
*/
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Used to STORE coefficients of Hermit interpolation polynoms */
/* KEYWORDS : */
/* ----------- */
/* HERMITE */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* The coefficients of hermit polynoms are calculated by */
/* the routine MMHERM0 and read by the routine MMHERM1 */
/* > */
/* **********************************************************************
*/
/* NBCOEF is the size of CMHERM (see below) */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* Data checking */
/* ***********************************************************************
*/
/* Parameter adjustments */
--debfin;
/* Function Body */
d1 = advapp_abs(debfin[1]);
if (d1 > (float)100.) {
goto L9101;
}
d2 = advapp_abs(debfin[2]);
if (d2 > (float)100.) {
goto L9101;
}
d2 = d1 + d2;
if (d2 < (float).01) {
goto L9101;
}
d1 = (d__1 = debfin[2] - debfin[1], advapp_abs(d__1));
if (d1 / d2 < (float).01) {
goto L9101;
}
/* ***********************************************************************
*/
/* Initialization */
/* ***********************************************************************
*/
*iercod = 0;
epspi = 1e-10;
/* ***********************************************************************
*/
/* IS IT ALREADY INITIALIZED ? */
d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
d1 *= 16111959;
if (debfin[1] != mmcmher_.tdebut) {
goto L100;
}
if (debfin[2] != mmcmher_.tfinal) {
goto L100;
}
if (d1 != mmcmher_.verifi) {
goto L100;
}
goto L9001;
/* ***********************************************************************
*/
/* CALCULATION */
/* ***********************************************************************
*/
L100:
/* Init. matrix identity : */
ncmat = 36;
AdvApp2Var_SysBase::mvriraz_(&ncmat,
miden);
for (ii = 1; ii <= 6; ++ii) {
miden[ii + ii * 6 - 7] = 1.;
/* L110: */
}
/* Init to 0 of table CMHERM */
AdvApp2Var_SysBase::mvriraz_(&c__576, mmcmher_.cmherm);
/* Calculation by solution of linear systems */
for (iord1 = -1; iord1 <= 2; ++iord1) {
for (iord2 = -1; iord2 <= 2; ++iord2) {
iord[0] = iord1;
iord[1] = iord2;
iof[0] = 0;
iof[1] = iord[0] + 1;
ncf = iord[0] + iord[1] + 2;
/* Calculate matrix MAT to invert: */
for (cot = 1; cot <= 2; ++cot) {
if (iord[cot - 1] > -1) {
prod = 1.;
i__1 = ncf;
for (jj = 1; jj <= i__1; ++jj) {
cof[jj - 1] = 1.;
/* L200: */
}
}
i__1 = iord[cot - 1] + 1;
for (pp = 1; pp <= i__1; ++pp) {
ii = pp + iof[cot - 1];
prod = 1.;
i__2 = pp - 1;
for (jj = 1; jj <= i__2; ++jj) {
mat[ii + jj * 6 - 7] = (float)0.;
/* L300: */
}
i__2 = ncf;
for (jj = pp; jj <= i__2; ++jj) {
/* everything is done in these 3 lines
*/
mat[ii + jj * 6 - 7] = cof[jj - 1] * prod;
cof[jj - 1] *= jj - pp;
prod *= debfin[cot];
/* L400: */
}
/* L500: */
}
/* L1000: */
}
/* Inversion */
if (ncf >= 1) {
AdvApp2Var_MathBase::mmmrslwd_(&c__6, &ncf, &ncf, mat, miden, &epspi, abid, amat, &
ier);
if (ier > 0) {
goto L9101;
}
}
for (cot = 1; cot <= 2; ++cot) {
i__1 = iord[cot - 1] + 1;
for (pp = 1; pp <= i__1; ++pp) {
i__2 = ncf;
for (ii = 1; ii <= i__2; ++ii) {
mmcmher_.cmherm[ii + (pp + (cot + ((iord1 + (iord2 <<
2)) << 1)) * 3) * 6 + 155] = amat[ii + (pp +
iof[cot - 1]) * 6 - 7];
/* L1300: */
}
/* L1400: */
}
/* L1500: */
}
/* L2000: */
}
/* L2010: */
}
/* ***********************************************************************
*/
/* The initialized flag is located: */
mmcmher_.tdebut = debfin[1];
mmcmher_.tfinal = debfin[2];
d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
mmcmher_.verifi = d1 * 16111959;
/* ***********************************************************************
*/
goto L9999;
/* ***********************************************************************
*/
L9101:
*iercod = 1;
goto L9999;
L9001:
*iercod = -1;
goto L9999;
/* ***********************************************************************
*/
L9999:
AdvApp2Var_SysBase::maermsg_("MMHERM0", iercod, 7L);
/* ***********************************************************************
*/
return 0 ;
} /* mmherm0_ */
//=======================================================================
//function : mmherm1_
//purpose :
//=======================================================================
int mmherm1_(doublereal *debfin,
integer *ordrmx,
integer *iordre,
doublereal *hermit,
integer *iercod)
{
/* System generated locals */
integer hermit_dim1, hermit_dim2, hermit_offset;
/* Local variables */
integer nbval;
doublereal d1;
integer cot;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* reading of coeffs. of HERMIT interpolation polynoms */
/* KEYWORDS : */
/* ----------- */
/* MATH_ACCES :: HERMIT */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* DEBFIN : PARAMETES DEFINING THE CONSTRAINTS */
/* DEBFIN(1) : FIRST PARAMETER */
/* DEBFIN(2) : SECOND PARAMETER */
/* Should be equal to the corresponding arguments during the */
/* last call to MMHERM0 for the initialization of coeffs. */
/* ORDRMX : indicates the dimensioning of HERMIT: */
/* there is no choice : ORDRMX should be equal to the value */
/* of PARAMETER IORDMX of INCLUDE MMCMHER, or 2 for the moment */
/* IORDRE (2) : Orders of constraints in each corresponding parameter DEBFIN(I) */
/* should be between -1 (no constraints) and ORDRMX. */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* HERMIT : HERMIT(1:IORDRE(1)+IORDRE(2)+2, j, cote) are the */
/* coefficients in the canonic base of Hermit polynom */
/* corresponding to orders IORDRE with parameters DEBFIN for */
/* the constraint of order j on DEBFIN(cote). j is between 0 and IORDRE(cote). */
/* IERCOD : Error code : */
/* -1: O.K but necessary to reinitialize the coefficients */
/* (info for optimization) */
/* 0 : O.K. */
/* 1 : Error in MMHERM0 */
/* 2 : arguments invalid */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Type Name */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* This program reads coefficients of Hermit polynoms */
/* that were earlier initialized by MMHERM0 */
/* PMN : initialisation is no more done by the caller. */
/* ***********************************************************************
*/
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Serves to STORE the coefficients of Hermit interpolation polynoms */
/* KEYWORDS : */
/* ----------- */
/* HERMITE */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* the coefficients of Hetmit polynoms are calculated by */
/* routine MMHERM0 and read by routine MMHERM1 */
/* > */
/* **********************************************************************
*/
/* NBCOEF is the size of CMHERM (see lower) */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* Initializations */
/* ***********************************************************************
*/
/* Parameter adjustments */
--debfin;
hermit_dim1 = (*ordrmx << 1) + 2;
hermit_dim2 = *ordrmx + 1;
hermit_offset = hermit_dim1 * hermit_dim2 + 1;
hermit -= hermit_offset;
--iordre;
/* Function Body */
*iercod = 0;
/* ***********************************************************************
*/
/* Data Checking */
/* ***********************************************************************
*/
if (*ordrmx != 2) {
goto L9102;
}
for (cot = 1; cot <= 2; ++cot) {
if (iordre[cot] < -1) {
goto L9102;
}
if (iordre[cot] > *ordrmx) {
goto L9102;
}
/* L100: */
}
/* IS-IT CORRECTLY INITIALIZED ? */
d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
d1 *= 16111959;
/* OTHERWISE IT IS INITIALIZED */
if (debfin[1] != mmcmher_.tdebut || debfin[2] != mmcmher_.tfinal || d1
!= mmcmher_.verifi) {
*iercod = -1;
mmherm0_(&debfin[1], iercod);
if (*iercod > 0) {
goto L9101;
}
}
/* ***********************************************************************
*/
/* READING */
/* ***********************************************************************
*/
nbval = 36;
AdvApp2Var_SysBase::msrfill_(&nbval, &mmcmher_.cmherm[((((iordre[1] + (iordre[2] << 2)) << 1)
+ 1) * 3 + 1) * 6 + 156], &hermit[hermit_offset]);
/* ***********************************************************************
*/
goto L9999;
/* ***********************************************************************
*/
L9101:
*iercod = 1;
goto L9999;
L9102:
*iercod = 2;
goto L9999;
/* ***********************************************************************
*/
L9999:
AdvApp2Var_SysBase::maermsg_("MMHERM1", iercod, 7L);
/* ***********************************************************************
*/
return 0 ;
} /* mmherm1_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmhjcan_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmhjcan_(integer *ndimen,
integer *ncourb,
integer *ncftab,
integer *orcont,
integer *ncflim,
doublereal *tcbold,
doublereal *tdecop,
doublereal *tcbnew,
integer *iercod)
{
integer c__2 = 2;
integer c__21 = 21;
/* System generated locals */
integer tcbold_dim1, tcbold_dim2, tcbold_offset, tcbnew_dim1, tcbnew_dim2,
tcbnew_offset, i__1, i__2, i__3, i__4, i__5;
/* Local variables */
logical ldbg;
integer ndeg;
doublereal taux1[21];
integer d__, e, i__, k;
doublereal mfact;
integer ncoeff;
doublereal tjacap[21];
integer iordre[2];
doublereal hermit[36]/* was [6][3][2] */, ctenor, bornes[2];
integer ier;
integer aux1, aux2;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* CONVERSION OF TABLE TCBOLD OF POLYNOMIAL CURVE COEFFICIENTS */
/* EXPRESSED IN HERMIT JACOBI BASE, INTO A */
/* TABLE OF COEFFICIENTS TCBNEW OF COURVES EXPRESSED IN THE CANONIC BASE */
/* KEYWORDS : */
/* ----------- */
/* CANNONIC, HERMIT, JACCOBI */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* ORDHER : ORDER OF HERMIT POLYNOMS OR ORDER OF CONTINUITY */
/* NCOEFS : NUMBER OF COEFFICIENTS OF A POLYNOMIAL CURVE */
/* FOR ONE OF ITS NDIM COMPONENTS;(DEGREE+1 OF THE CURVE)
*/
/* NDIM : DIMENSION OF THE CURVE */
/* CBHEJA : TABLE OF COEFFICIENTS OF THE CURVE IN THE BASE */
/* HERMIT JACOBI */
/* (H(0,-1),..,H(ORDHER,-1),H(0,1),..,H(ORDHER,1), */
/* JA(ORDHER+1,2*ORDHER+2),....,JA(ORDHER+1,NCOEFS-1) */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* CBRCAN : TABLE OF COEFFICIENTS OF THE CURVE IN THE CANONIC BASE */
/* (1, t, ...) */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Providesinteger constants from 0 to 1000 */
/* KEYWORDS : */
/* ----------- */
/* ALL, INTEGER */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* INITIALIZATION */
/* ***********************************************************************
*/
/* Parameter adjustments */
--ncftab;
tcbnew_dim1 = *ndimen;
tcbnew_dim2 = *ncflim;
tcbnew_offset = tcbnew_dim1 * (tcbnew_dim2 + 1) + 1;
tcbnew -= tcbnew_offset;
tcbold_dim1 = *ndimen;
tcbold_dim2 = *ncflim;
tcbold_offset = tcbold_dim1 * (tcbold_dim2 + 1) + 1;
tcbold -= tcbold_offset;
/* Function Body */
ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
if (ldbg) {
AdvApp2Var_SysBase::mgenmsg_("MMHJCAN", 7L);
}
*iercod = 0;
bornes[0] = -1.;
bornes[1] = 1.;
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
if (*orcont > 2) {
goto L9101;
}
if (*ncflim > 21) {
goto L9101;
}
/* CALCULATION OF HERMIT POLYNOMS IN THE CANONIC BASE ON (-1,1) */
iordre[0] = *orcont;
iordre[1] = *orcont;
mmherm1_(bornes, &c__2, iordre, hermit, &ier);
if (ier > 0) {
goto L9102;
}
aux1 = *orcont + 1;
aux2 = aux1 << 1;
i__1 = *ncourb;
for (e = 1; e <= i__1; ++e) {
ctenor = (tdecop[e] - tdecop[e - 1]) / 2;
ncoeff = ncftab[e];
ndeg = ncoeff - 1;
if (ncoeff > 21) {
goto L9101;
}
i__2 = *ndimen;
for (d__ = 1; d__ <= i__2; ++d__) {
/* CONVERSION OF THE COEFFICIENTS OF THE PART OF THE CURVE EXPRESSED */
/* IN HERMIT BASE, INTO THE CANONIC BASE */
AdvApp2Var_SysBase::mvriraz_(&ncoeff, taux1);
i__3 = aux2;
for (k = 1; k <= i__3; ++k) {
i__4 = aux1;
for (i__ = 1; i__ <= i__4; ++i__) {
i__5 = i__ - 1;
mfact = AdvApp2Var_MathBase::pow__di(&ctenor, &i__5);
taux1[k - 1] += (tcbold[d__ + (i__ + e * tcbold_dim2) *
tcbold_dim1] * hermit[k + (i__ + 2) * 6 - 19] +
tcbold[d__ + (i__ + aux1 + e * tcbold_dim2) *
tcbold_dim1] * hermit[k + (i__ + 5) * 6 - 19]) *
mfact;
}
}
i__3 = ncoeff;
for (i__ = aux2 + 1; i__ <= i__3; ++i__) {
taux1[i__ - 1] = tcbold[d__ + (i__ + e * tcbold_dim2) *
tcbold_dim1];
}
/* CONVERSION OF THE COEFFICIENTS OF THE PART OF THE CURVE EXPRESSED */
/* IN CANONIC-JACOBI BASE, INTO THE CANONIC BASE */
AdvApp2Var_MathBase::mmapcmp_(&minombr_.nbr[1], &c__21, &ncoeff, taux1, tjacap);
AdvApp2Var_MathBase::mmjacan_(orcont, &ndeg, tjacap, taux1);
/* RECOPY THE COEFS RESULTING FROM THE CONVERSION IN THE TABLE */
/* OF RESULTS */
i__3 = ncoeff;
for (i__ = 1; i__ <= i__3; ++i__) {
tcbnew[d__ + (i__ + e * tcbnew_dim2) * tcbnew_dim1] = taux1[
i__ - 1];
}
}
}
goto L9999;
/* ***********************************************************************
*/
/* PROCESSING OF ERRORS */
/* ***********************************************************************
*/
L9101:
*iercod = 1;
goto L9999;
L9102:
*iercod = 2;
goto L9999;
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
L9999:
AdvApp2Var_SysBase::maermsg_("MMHJCAN", iercod, 7L);
if (ldbg) {
AdvApp2Var_SysBase::mgsomsg_("MMHJCAN", 7L);
}
return 0 ;
} /* mmhjcan_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mminltt_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mminltt_(integer *ncolmx,
integer *nlgnmx,
doublereal *tabtri,
integer *nbrcol,
integer *nbrlgn,
doublereal *ajoute,
doublereal *,//epseg,
integer *iercod)
{
/* System generated locals */
integer tabtri_dim1, tabtri_offset, i__1, i__2;
/* Local variables */
logical idbg;
integer icol, ilgn, nlgn, noct, inser;
doublereal epsega = 0.;
integer ibb;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* . Insert a line in a table parsed without redundance */
/* KEYWORDS : */
/* ----------- */
/* TOUS,MATH_ACCES :: TABLEAU&,INSERTION,&TABLEAU */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* . NCOLMX : Number of columns in the table */
/* . NLGNMX : Number of lines in the table */
/* . TABTRI : Table parsed by lines without redundances */
/* . NBRCOL : Number of columns used */
/* . NBRLGN : Number of lines used */
/* . AJOUTE : Line to be added */
/* . EPSEGA : Epsilon to test the redundance */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* . TABTRI : Table parsed by lines without redundances */
/* . NBRLGN : Number of lines used */
/* . IERCOD : 0 -> No problem */
/* 1 -> The table is full */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* . The line is inserted only if there is no line with all
*/
/* elements equl to those which are planned to be insered, to epsilon. */
/* . Level of de debug = 3 */
/**/
/* DECLARATIONS , CONTROL OF INPUT ARGUMENTS , INITIALIZATION */
/* ***********************************************************************
*/
/* --- Parameters */
/* --- Functions */
/* --- Local variables */
/* --- Messages */
/* Parameter adjustments */
tabtri_dim1 = *ncolmx;
tabtri_offset = tabtri_dim1 + 1;
tabtri -= tabtri_offset;
--ajoute;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
idbg = ibb >= 3;
if (idbg) {
AdvApp2Var_SysBase::mgenmsg_("MMINLTT", 7L);
}
/* --- Control arguments */
if (*nbrlgn >= *nlgnmx) {
goto L9001;
}
/* -------------------- */
/* *** INITIALIZATION */
/* -------------------- */
*iercod = 0;
/* ---------------------------- */
/* *** SEARCH OF REDUNDANCE */
/* ---------------------------- */
i__1 = *nbrlgn;
for (ilgn = 1; ilgn <= i__1; ++ilgn) {
if (tabtri[ilgn * tabtri_dim1 + 1] >= ajoute[1] - epsega) {
if (tabtri[ilgn * tabtri_dim1 + 1] <= ajoute[1] + epsega) {
i__2 = *nbrcol;
for (icol = 1; icol <= i__2; ++icol) {
if (tabtri[icol + ilgn * tabtri_dim1] < ajoute[icol] -
epsega || tabtri[icol + ilgn * tabtri_dim1] >
ajoute[icol] + epsega) {
goto L20;
}
/* L10: */
}
goto L9999;
} else {
goto L30;
}
}
L20:
;
}
/* ----------------------------------- */
/* *** SEARCH OF THE INSERTION POINT */
/* ----------------------------------- */
L30:
i__1 = *nbrlgn;
for (ilgn = 1; ilgn <= i__1; ++ilgn) {
i__2 = *nbrcol;
for (icol = 1; icol <= i__2; ++icol) {
if (tabtri[icol + ilgn * tabtri_dim1] < ajoute[icol]) {
goto L50;
}
if (tabtri[icol + ilgn * tabtri_dim1] > ajoute[icol]) {
goto L70;
}
/* L60: */
}
L50:
;
}
ilgn = *nbrlgn + 1;
/* -------------- */
/* *** INSERTION */
/* -------------- */
L70:
inser = ilgn;
++(*nbrlgn);
/* --- Shift lower */
nlgn = *nbrlgn - inser;
if (nlgn > 0) {
noct = (*ncolmx << 3) * nlgn;
AdvApp2Var_SysBase::mcrfill_(&noct,
&tabtri[inser * tabtri_dim1 + 1],
&tabtri[(inser + 1)* tabtri_dim1 + 1]);
}
/* --- Copy line */
noct = *nbrcol << 3;
AdvApp2Var_SysBase::mcrfill_(&noct,
&ajoute[1],
&tabtri[inser * tabtri_dim1 + 1]);
goto L9999;
/* ******************************************************************** */
/* OUTPUT ERROR , RETURN CALLING PROGRAM , MESSAGES */
/* ******************************************************************** */
/* --- The table is already full */
L9001:
*iercod = 1;
/* --- End */
L9999:
if (*iercod != 0) {
AdvApp2Var_SysBase::maermsg_("MMINLTT", iercod, 7L);
}
if (idbg) {
AdvApp2Var_SysBase::mgsomsg_("MMINLTT", 7L);
}
return 0 ;
} /* mminltt_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmjacan_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmjacan_(const integer *ideriv,
integer *ndeg,
doublereal *poljac,
doublereal *polcan)
{
/* System generated locals */
integer poljac_dim1, i__1, i__2;
/* Local variables */
integer iptt, i__, j, ibb;
doublereal bid;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Routine of transfer of Jacobi normalized to canonic [-1,1], */
/* the tables are ranked by even, then by uneven degree. */
/* KEYWORDS : */
/* ----------- */
/* LEGENDRE,JACOBI,PASSAGE. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* IDERIV : Order of Jacobi between -1 and 2. */
/* NDEG : The true degree of the polynom. */
/* POLJAC : The polynom in the Jacobi base. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* POLCAN : The curve expressed in the canonic base [-1,1]. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* Name of the routine */
/* Matrices of conversion */
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* MATRIX OF TRANSFORMATION OF LEGENDRE BASE */
/* KEYWORDS : */
/* ----------- */
/* MATH */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* Legendre common / Restricted Casteljau. */
/* 0:1 0 Concerns the even terms, 1 the uneven terms. */
/* CANPLG : Matrix of passage to canonic from Jacobi with calculated parities */
/* PLGCAN : Matrix of passage from Jacobi to canonic with calculated parities */
/* ***********************************************************************
*/
/* Parameter adjustments */
poljac_dim1 = *ndeg / 2 + 1;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 5) {
AdvApp2Var_SysBase::mgenmsg_("MMJACAN", 7L);
}
/* ----------------- Expression of terms of even degree ----------------
*/
i__1 = *ndeg / 2;
for (i__ = 0; i__ <= i__1; ++i__) {
bid = 0.;
iptt = i__ * 31 - (i__ + 1) * i__ / 2 + 1;
i__2 = *ndeg / 2;
for (j = i__; j <= i__2; ++j) {
bid += mmjcobi_.plgcan[iptt + j + *ideriv * 992 + 991] * poljac[
j];
/* L310: */
}
polcan[i__ * 2] = bid;
/* L300: */
}
/* --------------- Expression of terms of uneven degree ----------------
*/
if (*ndeg == 0) {
goto L9999;
}
i__1 = (*ndeg - 1) / 2;
for (i__ = 0; i__ <= i__1; ++i__) {
bid = 0.;
iptt = i__ * 31 - (i__ + 1) * i__ / 2 + 1;
i__2 = (*ndeg - 1) / 2;
for (j = i__; j <= i__2; ++j) {
bid += mmjcobi_.plgcan[iptt + j + ((*ideriv << 1) + 1) * 496 +
991] * poljac[j + poljac_dim1];
/* L410: */
}
polcan[(i__ << 1) + 1] = bid;
/* L400: */
}
/* -------------------------------- The end -----------------------------
*/
L9999:
if (ibb >= 5) {
AdvApp2Var_SysBase::mgsomsg_("MMJACAN", 7L);
}
return 0;
} /* mmjacan_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmjaccv_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmjaccv_(const integer *ncoef,
const integer *ndim,
const integer *ider,
const doublereal *crvlgd,
doublereal *polaux,
doublereal *crvcan)
{
/* Initialized data */
static char nomprg[8+1] = "MMJACCV ";
/* System generated locals */
integer crvlgd_dim1, crvlgd_offset, crvcan_dim1, crvcan_offset,
polaux_dim1, i__1, i__2;
/* Local variables */
integer ndeg, i__, nd, ii, ibb;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Passage from the normalized Jacobi base to the canonic base. */
/* KEYWORDS : */
/* ----------- */
/* SMOOTHING, BASE, LEGENDRE */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIM: Space Dimension. */
/* NCOEF: Degree +1 of the polynom. */
/* IDER: Order of Jacobi polynoms. */
/* CRVLGD : Curve in the base of Jacobi. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* POLAUX : Auxilliary space. */
/* CRVCAN : The curve in the canonic base [-1,1] */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* *********************************************************************
*/
/* Name of the routine */
/* Parameter adjustments */
polaux_dim1 = (*ncoef - 1) / 2 + 1;
crvcan_dim1 = *ncoef - 1 + 1;
crvcan_offset = crvcan_dim1;
crvcan -= crvcan_offset;
crvlgd_dim1 = *ncoef - 1 + 1;
crvlgd_offset = crvlgd_dim1;
crvlgd -= crvlgd_offset;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 3) {
AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
}
ndeg = *ncoef - 1;
i__1 = *ndim;
for (nd = 1; nd <= i__1; ++nd) {
/* Loading of the auxilliary table. */
ii = 0;
i__2 = ndeg / 2;
for (i__ = 0; i__ <= i__2; ++i__) {
polaux[i__] = crvlgd[ii + nd * crvlgd_dim1];
ii += 2;
/* L310: */
}
ii = 1;
if (ndeg >= 1) {
i__2 = (ndeg - 1) / 2;
for (i__ = 0; i__ <= i__2; ++i__) {
polaux[i__ + polaux_dim1] = crvlgd[ii + nd * crvlgd_dim1];
ii += 2;
/* L320: */
}
}
/* Call the routine of base change. */
AdvApp2Var_MathBase::mmjacan_(ider, &ndeg, polaux, &crvcan[nd * crvcan_dim1]);
/* L300: */
}
/* L9999: */
return 0;
} /* mmjaccv_ */
//=======================================================================
//function : mmloncv_
//purpose :
//=======================================================================
int mmloncv_(integer *ndimax,
integer *ndimen,
integer *ncoeff,
doublereal *courbe,
doublereal *tdebut,
doublereal *tfinal,
doublereal *xlongc,
integer *iercod)
{
/* Initialized data */
integer kgar = 0;
/* System generated locals */
integer courbe_dim1, courbe_offset, i__1, i__2;
/* Local variables */
doublereal tran;
integer ngaus = 0;
doublereal c1, c2, d1, d2,
wgaus[20] = {0.}, uroot[20] = {0.}, x1, x2, dd;
integer ii, jj, kk;
doublereal som;
doublereal der1, der2;
/* **********************************************************************
*/
/* FUNCTION : Length of an arc of curve on a given interval */
/* ---------- for a function the mathematic representation */
/* which of is a multidimensional polynom. */
/* The polynom is a set of polynoms the coefficients which of are ranked */
/* in a table with 2 indices, each line relative to 1 polynom. */
/* The polynom is defined by its coefficients ordered by increasing
* power of the variable. */
/* All polynoms have the same number of coefficients (and the same degree). */
/* KEYWORDS : LENGTH, CURVE */
/* ----------- */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* NDIMAX : Max number of lines of tables (max number of polynoms). */
/* NDIMEN : Dimension of the polynom (Nomber of polynoms). */
/* NCOEFF : Number of coefficients of the polynom (no limitation) */
/* This is degree + 1 */
/* COURBE : Coefficients of the polynom ordered by increasing power */
/* Dimension to (NDIMAX,NCOEFF). */
/* TDEBUT : Lower limit of integration for length calculation. */
/* TFINAL : Upper limit of integration for length calculation. */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* XLONGC : Length of arc of curve */
/* IERCOD : Error code : */
/* = 0 ==> All is OK */
/* = 1 ==> NDIMEN or NCOEFF negative or null */
/* = 2 ==> Pb loading Legendre roots and Gauss weight */
/* by MVGAUS0. */
/* If error => XLONGC = 0 */
/* COMMONS USED : */
/* ------------------ */
/* .Neant. */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Type Name */
/* MAERMSG R*8 DSQRT I*4 MIN */
/* MVGAUS0 */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* See VGAUSS to understand well the technique. */
/* Actually SQRT (dpi^2) is integrated for i=1,nbdime */
/* Calculation of the derivative is included in the code to avoid an additional */
/* call of the routine. */
/* The integrated function is strictly increasing, it */
/* is not necessary to use a high degree for the GAUSS method GAUSS. */
/* The degree of LEGENDRE polynom results from the degree of the */
/* polynom to be integrated. It can vary from 4 to 40 (with step of 4). */
/* The precision (relative) of integration is of order 1.D-8. */
/* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
/* Attention : the precision of the result is not controlled. */
/* If you wish to control it, use MMCGLC1, taking into account that */
/* the performance (in time) will be worse. */
/* >=====================================================================
*/
/* ATTENTION : SAVE KGAR WGAUS and UROOT EVENTUALLY */
/* ,IERXV */
/* INTEGER I1,I20 */
/* PARAMETER (I1=1,I20=20) */
/* Parameter adjustments */
courbe_dim1 = *ndimax;
courbe_offset = courbe_dim1 + 1;
courbe -= courbe_offset;
/* Function Body */
/* ****** General initialization ** */
*iercod = 999999;
*xlongc = 0.;
/* ****** Initialization of UROOT, WGAUS, NGAUS and KGAR ** */
/* CALL MXVINIT(IERXV,'INTEGER',I1,KGAR,'INTEGER',I1,NGAUS */
/* 1 ,'DOUBLE PRECISION',I20,UROOT,'DOUBLE PRECISION',I20,WGAUS) */
/* IF (IERXV.GT.0) KGAR=0 */
/* ****** Test the equity of limits ** */
if (*tdebut == *tfinal) {
*iercod = 0;
goto L9900;
}
/* ****** Test the dimension and the number of coefficients ** */
if (*ndimen <= 0 || *ncoeff <= 0) {
*iercod = 1;
goto L9900;
}
/* ****** Calculate the optimal degree ** */
kk = *ncoeff / 4 + 1;
kk = advapp_min(kk,10);
/* ****** Return the coefficients for the integral (DEGRE=4*KK) */
/* if KK <> KGAR. */
if (kk != kgar) {
mvgaus0_(&kk, uroot, wgaus, &ngaus, iercod);
if (*iercod > 0) {
kgar = 0;
*iercod = 2;
goto L9900;
}
kgar = kk;
}
/* C1 => Point medium interval */
/* C2 => 1/2 amplitude interval */
c1 = (*tfinal + *tdebut) * .5;
c2 = (*tfinal - *tdebut) * .5;
/* ----------------------------------------------------------- */
/* ****** Integration - Loop on GAUSS intervals ** */
/* ----------------------------------------------------------- */
som = 0.;
i__1 = ngaus;
for (jj = 1; jj <= i__1; ++jj) {
/* ****** Integration taking the symmetry into account ** */
tran = c2 * uroot[jj - 1];
x1 = c1 + tran;
x2 = c1 - tran;
/* ****** Derivation on the dimension of the space ** */
der1 = 0.;
der2 = 0.;
i__2 = *ndimen;
for (kk = 1; kk <= i__2; ++kk) {
d1 = (*ncoeff - 1) * courbe[kk + *ncoeff * courbe_dim1];
d2 = d1;
for (ii = *ncoeff - 1; ii >= 2; --ii) {
dd = (ii - 1) * courbe[kk + ii * courbe_dim1];
d1 = d1 * x1 + dd;
d2 = d2 * x2 + dd;
/* L100: */
}
der1 += d1 * d1;
der2 += d2 * d2;
/* L200: */
}
/* ****** Integration ** */
som += wgaus[jj - 1] * c2 * (sqrt(der1) + sqrt(der2));
/* ****** End of loop on GAUSS intervals ** */
/* L300: */
}
/* ****** Work ended ** */
*xlongc = som;
/* ****** It is forced IERCOD = 0 ** */
*iercod = 0;
/* ****** Final processing ** */
L9900:
/* ****** Save UROOT, WGAUS, NGAUS and KGAR ** */
/* CALL MXVSAVE(IERXV,'INTEGER',I1,KGAR,'INTEGER',I1,NGAUS */
/* 1 ,'DOUBLE PRECISION',I20,UROOT,'DOUBLE PRECISION',I20,WGAUS) */
/* IF (IERXV.GT.0) KGAR=0 */
/* ****** End of sub-program ** */
if (*iercod != 0) {
AdvApp2Var_SysBase::maermsg_("MMLONCV", iercod, 7L);
}
return 0 ;
} /* mmloncv_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmpobas_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmpobas_(doublereal *tparam,
integer *iordre,
integer *ncoeff,
integer *nderiv,
doublereal *valbas,
integer *iercod)
{
integer c__2 = 2;
integer c__1 = 1;
/* Initialized data */
doublereal moin11[2] = { -1.,1. };
/* System generated locals */
integer valbas_dim1, i__1;
/* Local variables */
doublereal vjac[80], herm[24];
integer iord[2];
doublereal wval[4];
integer nwcof, iunit;
doublereal wpoly[7];
integer ii, jj, iorjac;
doublereal hermit[36] /* was [6][3][2] */;
integer kk1, kk2, kk3;
integer khe, ier;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Position on the polynoms of base hermit-Jacobi */
/* and their succesive derivatives */
/* KEYWORDS : */
/* ----------- */
/* PUBLIC, POSITION, HERMIT, JACOBI */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* TPARAM : Parameter for which the position is found. */
/* IORDRE : Orderof hermit-Jacobi (-1,0,1, ou 2) */
/* NCOEFF : Number of coefficients of polynoms (Nb of value to calculate) */
/* NDERIV : Number of derivative to calculate (0<= N <=3) */
/* 0 -> Position simple on base functions */
/* N -> Position on base functions and derivative */
/* of order 1 to N */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* VALBAS (NCOEFF, 0:NDERIV) : calculated value */
/* i */
/* d vj(t) = VALBAS(J, I) */
/* -- i */
/* dt */
/* IERCOD : Error code */
/* 0 : Ok */
/* 1 : Incoherence of input arguments */
/* COMMONS USED : */
/* -------------- */
/* REFERENCES CALLED : */
/* ------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* Parameter adjustments */
valbas_dim1 = *ncoeff;
--valbas;
/* Function Body */
/* ***********************************************************************
*/
/* INITIALIZATIONS */
/* ***********************************************************************
*/
*iercod = 0;
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
if (*nderiv > 3) {
goto L9101;
}
if (*ncoeff > 20) {
goto L9101;
}
if (*iordre > 2) {
goto L9101;
}
iord[0] = *iordre;
iord[1] = *iordre;
iorjac = (*iordre + 1) << 1;
/* (1) Generic Calculations .... */
/* (1.a) Calculation of hermit polynoms */
if (*iordre >= 0) {
mmherm1_(moin11, &c__2, iord, hermit, &ier);
if (ier > 0) {
goto L9102;
}
}
/* (1.b) Evaluation of hermit polynoms */
jj = 1;
iunit = *nderiv + 1;
khe = (*iordre + 1) * iunit;
if (*nderiv > 0) {
i__1 = *iordre;
for (ii = 0; ii <= i__1; ++ii) {
mmdrvcb_(nderiv, &c__1, &iorjac, &hermit[(ii + 3) * 6 - 18],
tparam, &herm[jj - 1], &ier);
if (ier > 0) {
goto L9102;
}
mmdrvcb_(nderiv, &c__1, &iorjac, &hermit[(ii + 6) * 6 - 18],
tparam, &herm[jj + khe - 1], &ier);
if (ier > 0) {
goto L9102;
}
jj += iunit;
}
} else {
i__1 = *iordre;
for (ii = 0; ii <= i__1; ++ii) {
AdvApp2Var_MathBase::mmpocrb_(&c__1, &iorjac, &hermit[(ii + 3) * 6 - 18], &c__1,
tparam, &herm[jj - 1]);
AdvApp2Var_MathBase::mmpocrb_(&c__1, &iorjac, &hermit[(ii + 6) * 6 - 18], &c__1,
tparam, &herm[jj + khe - 1]);
jj += iunit;
}
}
/* (1.c) Evaluation of Jacobi polynoms */
ii = *ncoeff - iorjac;
mmpojac_(tparam, &iorjac, &ii, nderiv, vjac, &ier);
if (ier > 0) {
goto L9102;
}
/* (1.d) Evaluation of W(t) */
/* Computing MAX */
i__1 = iorjac + 1;
nwcof = advapp_max(i__1,1);
AdvApp2Var_SysBase::mvriraz_(&nwcof,
wpoly);
wpoly[0] = 1.;
if (*iordre == 2) {
wpoly[2] = -3.;
wpoly[4] = 3.;
wpoly[6] = -1.;
} else if (*iordre == 1) {
wpoly[2] = -2.;
wpoly[4] = 1.;
} else if (*iordre == 0) {
wpoly[2] = -1.;
}
mmdrvcb_(nderiv, &c__1, &nwcof, wpoly, tparam, wval, &ier);
if (ier > 0) {
goto L9102;
}
kk1 = *ncoeff - iorjac;
kk2 = kk1 << 1;
kk3 = kk1 * 3;
/* (2) Evaluation of order 0 */
jj = 1;
i__1 = iorjac;
for (ii = 1; ii <= i__1; ++ii) {
valbas[ii] = herm[jj - 1];
jj += iunit;
}
i__1 = kk1;
for (ii = 1; ii <= i__1; ++ii) {
valbas[ii + iorjac] = wval[0] * vjac[ii - 1];
}
/* (3) Evaluation of order 1 */
if (*nderiv >= 1) {
jj = 2;
i__1 = iorjac;
for (ii = 1; ii <= i__1; ++ii) {
valbas[ii + valbas_dim1] = herm[jj - 1];
jj += iunit;
}
i__1 = kk1;
for (ii = 1; ii <= i__1; ++ii) {
valbas[ii + iorjac + valbas_dim1] = wval[0] * vjac[ii + kk1 - 1]
+ wval[1] * vjac[ii - 1];
}
}
/* (4) Evaluation of order 2 */
if (*nderiv >= 2) {
jj = 3;
i__1 = iorjac;
for (ii = 1; ii <= i__1; ++ii) {
valbas[ii + (valbas_dim1 << 1)] = herm[jj - 1];
jj += iunit;
}
i__1 = kk1;
for (ii = 1; ii <= i__1; ++ii) {
valbas[ii + iorjac + (valbas_dim1 << 1)] = wval[0] * vjac[ii +
kk2 - 1] + wval[1] * 2 * vjac[ii + kk1 - 1] + wval[2] *
vjac[ii - 1];
}
}
/* (5) Evaluation of order 3 */
if (*nderiv >= 3) {
jj = 4;
i__1 = iorjac;
for (ii = 1; ii <= i__1; ++ii) {
valbas[ii + valbas_dim1 * 3] = herm[jj - 1];
jj += iunit;
}
i__1 = kk1;
for (ii = 1; ii <= i__1; ++ii) {
valbas[ii + iorjac + valbas_dim1 * 3] = wval[0] * vjac[ii + kk3 -
1] + wval[1] * 3 * vjac[ii + kk2 - 1] + wval[2] * 3 *
vjac[ii + kk1 - 1] + wval[3] * vjac[ii - 1];
}
}
goto L9999;
/* ***********************************************************************
*/
/* ERROR PROCESSING */
/* ***********************************************************************
*/
L9101:
*iercod = 1;
goto L9999;
L9102:
*iercod = 2;
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
L9999:
if (*iercod > 0) {
AdvApp2Var_SysBase::maermsg_("MMPOBAS", iercod, 7L);
}
return 0 ;
} /* mmpobas_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmpocrb_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmpocrb_(integer *ndimax,
integer *ncoeff,
doublereal *courbe,
integer *ndim,
doublereal *tparam,
doublereal *pntcrb)
{
/* System generated locals */
integer courbe_dim1, courbe_offset, i__1, i__2;
/* Local variables */
integer ncof2;
integer isize, nd, kcf, ncf;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* CALCULATE THE COORDINATES OF A POINT OF A CURVE OF GIVEN PARAMETER */
/* TPARAM ( IN 2D, 3D OR MORE) */
/* KEYWORDS : */
/* ----------- */
/* TOUS , MATH_ACCES :: COURBE&,PARAMETRE& , POSITIONNEMENT , &POINT
*/
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMAX : format / dimension of the curve */
/* NCOEFF : Nb of coefficients of the curve */
/* COURBE : Matrix of coefficients of the curve */
/* NDIM : Dimension useful of the workspace */
/* TPARAM : Value of the parameter where the point is calculated */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* PNTCRB : Coordinates of the calculated point */
/* COMMONS USED : */
/* ---------------- */
/* .Neant. */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Type Name */
/* MIRAZ MVPSCR2 MVPSCR3 */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
courbe_dim1 = *ndimax;
courbe_offset = courbe_dim1 + 1;
courbe -= courbe_offset;
--pntcrb;
/* Function Body */
isize = *ndim << 3;
AdvApp2Var_SysBase::miraz_(&isize,
&pntcrb[1]);
if (*ncoeff <= 0) {
goto L9999;
}
/* optimal processing 3d */
if (*ndim == 3 && *ndimax == 3) {
mvpscr3_(ncoeff, &courbe[courbe_offset], tparam, &pntcrb[1]);
/* optimal processing 2d */
} else if (*ndim == 2 && *ndimax == 2) {
mvpscr2_(ncoeff, &courbe[courbe_offset], tparam, &pntcrb[1]);
/* Any dimension - scheme of HORNER */
} else if (*tparam == 0.) {
i__1 = *ndim;
for (nd = 1; nd <= i__1; ++nd) {
pntcrb[nd] = courbe[nd + courbe_dim1];
/* L100: */
}
} else if (*tparam == 1.) {
i__1 = *ncoeff;
for (ncf = 1; ncf <= i__1; ++ncf) {
i__2 = *ndim;
for (nd = 1; nd <= i__2; ++nd) {
pntcrb[nd] += courbe[nd + ncf * courbe_dim1];
/* L300: */
}
/* L200: */
}
} else {
ncof2 = *ncoeff + 2;
i__1 = *ndim;
for (nd = 1; nd <= i__1; ++nd) {
i__2 = *ncoeff;
for (ncf = 2; ncf <= i__2; ++ncf) {
kcf = ncof2 - ncf;
pntcrb[nd] = (pntcrb[nd] + courbe[nd + kcf * courbe_dim1]) * *
tparam;
/* L500: */
}
pntcrb[nd] += courbe[nd + courbe_dim1];
/* L400: */
}
}
L9999:
return 0 ;
} /* mmpocrb_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmmpocur_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmmpocur_(integer *ncofmx,
integer *ndim,
integer *ndeg,
doublereal *courbe,
doublereal *tparam,
doublereal *tabval)
{
/* System generated locals */
integer courbe_dim1, courbe_offset, i__1;
/* Local variables */
integer i__, nd;
doublereal fu;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Position of a point on curve (ncofmx,ndim). */
/* KEYWORDS : */
/* ----------- */
/* TOUS , AB_SPECIFI :: COURBE&,POLYNOME&,POSITIONNEMENT,&POINT */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX: Format / degree of the CURVE. */
/* NDIM : Dimension of the space. */
/* NDEG : Degree of the polynom. */
/* COURBE: Coefficients of the curve. */
/* TPARAM: Parameter on the curve */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* TABVAL(NDIM): The resulting point (or table of values) */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* Parameter adjustments */
--tabval;
courbe_dim1 = *ncofmx;
courbe_offset = courbe_dim1 + 1;
courbe -= courbe_offset;
/* Function Body */
if (*ndeg < 1) {
i__1 = *ndim;
for (nd = 1; nd <= i__1; ++nd) {
tabval[nd] = 0.;
/* L290: */
}
} else {
i__1 = *ndim;
for (nd = 1; nd <= i__1; ++nd) {
fu = courbe[*ndeg + nd * courbe_dim1];
for (i__ = *ndeg - 1; i__ >= 1; --i__) {
fu = fu * *tparam + courbe[i__ + nd * courbe_dim1];
/* L120: */
}
tabval[nd] = fu;
/* L300: */
}
}
return 0 ;
} /* mmmpocur_ */
//=======================================================================
//function : mmpojac_
//purpose :
//=======================================================================
int mmpojac_(doublereal *tparam,
integer *iordre,
integer *ncoeff,
integer *nderiv,
doublereal *valjac,
integer *iercod)
{
integer c__2 = 2;
/* System generated locals */
integer valjac_dim1, i__1, i__2;
/* Local variables */
doublereal cofa, cofb, denom, tnorm[100];
integer ii, jj, kk1, kk2;
doublereal aux1, aux2;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Positioning on Jacobi polynoms and their derivatives */
/* successive by a recurrent algorithm */
/* KEYWORDS : */
/* ----------- */
/* RESERVE, POSITIONING, JACOBI */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* TPARAM : Parameter for which positioning is done. */
/* IORDRE : Order of hermit-?? (-1,0,1, or 2) */
/* NCOEFF : Number of coeeficients of polynoms (Nb of value to */
/* calculate) */
/* NDERIV : Number of derivative to calculate (0<= N <=3) */
/* 0 -> Position simple on jacobi functions */
/* N -> Position on jacobi functions and their */
/* derivatives of order 1 to N. */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* VALJAC (NCOEFF, 0:NDERIV) : the calculated values */
/* i */
/* d vj(t) = VALJAC(J, I) */
/* -- i */
/* dt */
/* IERCOD : Error Code */
/* 0 : Ok */
/* 1 : Incoherence of input arguments */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* static varaibles */
/* Parameter adjustments */
valjac_dim1 = *ncoeff;
--valjac;
/* Function Body */
/* ***********************************************************************
*/
/* INITIALISATIONS */
/* ***********************************************************************
*/
*iercod = 0;
/* ***********************************************************************
*/
/* Processing */
/* ***********************************************************************
*/
if (*nderiv > 3) {
goto L9101;
}
if (*ncoeff > 100) {
goto L9101;
}
/* --- Calculation of norms */
/* IF (NCOEFF.GT.NBCOF) THEN */
i__1 = *ncoeff;
for (ii = 1; ii <= i__1; ++ii) {
kk1 = ii - 1;
aux2 = 1.;
i__2 = *iordre;
for (jj = 1; jj <= i__2; ++jj) {
aux2 = aux2 * (doublereal) (kk1 + *iordre + jj) / (doublereal) (
kk1 + jj);
}
i__2 = (*iordre << 1) + 1;
tnorm[ii - 1] = sqrt(aux2 * (kk1 * 2. + (*iordre << 1) + 1) / pow__ii(&
c__2, &i__2));
}
/* END IF */
/* --- Trivial Positions ----- */
valjac[1] = 1.;
aux1 = (doublereal) (*iordre + 1);
valjac[2] = aux1 * *tparam;
if (*nderiv >= 1) {
valjac[valjac_dim1 + 1] = 0.;
valjac[valjac_dim1 + 2] = aux1;
if (*nderiv >= 2) {
valjac[(valjac_dim1 << 1) + 1] = 0.;
valjac[(valjac_dim1 << 1) + 2] = 0.;
if (*nderiv >= 3) {
valjac[valjac_dim1 * 3 + 1] = 0.;
valjac[valjac_dim1 * 3 + 2] = 0.;
}
}
}
/* --- Positioning by recurrence */
i__1 = *ncoeff;
for (ii = 3; ii <= i__1; ++ii) {
kk1 = ii - 1;
kk2 = ii - 2;
aux1 = (doublereal) (*iordre + kk2);
aux2 = aux1 * 2;
cofa = aux2 * (aux2 + 1) * (aux2 + 2);
cofb = (aux2 + 2) * -2. * aux1 * aux1;
denom = kk1 * 2. * (kk2 + (*iordre << 1) + 1) * aux2;
denom = 1. / denom;
/* --> Pi(t) */
valjac[ii] = (cofa * *tparam * valjac[kk1] + cofb * valjac[kk2]) *
denom;
/* --> P'i(t) */
if (*nderiv >= 1) {
valjac[ii + valjac_dim1] = (cofa * *tparam * valjac[kk1 +
valjac_dim1] + cofa * valjac[kk1] + cofb * valjac[kk2 +
valjac_dim1]) * denom;
/* --> P''i(t) */
if (*nderiv >= 2) {
valjac[ii + (valjac_dim1 << 1)] = (cofa * *tparam * valjac[
kk1 + (valjac_dim1 << 1)] + cofa * 2 * valjac[kk1 +
valjac_dim1] + cofb * valjac[kk2 + (valjac_dim1 << 1)]
) * denom;
}
/* --> P'i(t) */
if (*nderiv >= 3) {
valjac[ii + valjac_dim1 * 3] = (cofa * *tparam * valjac[kk1 +
valjac_dim1 * 3] + cofa * 3 * valjac[kk1 + (
valjac_dim1 << 1)] + cofb * valjac[kk2 + valjac_dim1 *
3]) * denom;
}
}
}
/* ---> Normalization */
i__1 = *ncoeff;
for (ii = 1; ii <= i__1; ++ii) {
i__2 = *nderiv;
for (jj = 0; jj <= i__2; ++jj) {
valjac[ii + jj * valjac_dim1] = tnorm[ii - 1] * valjac[ii + jj *
valjac_dim1];
}
}
goto L9999;
/* ***********************************************************************
*/
/* PROCESSING OF ERRORS */
/* ***********************************************************************
*/
L9101:
*iercod = 1;
goto L9999;
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
L9999:
if (*iercod > 0) {
AdvApp2Var_SysBase::maermsg_("MMPOJAC", iercod, 7L);
}
return 0 ;
} /* mmpojac_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmposui_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmposui_(integer *dimmat,
integer *,//nistoc,
integer *aposit,
integer *posuiv,
integer *iercod)
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
logical ldbg;
integer imin, jmin, i__, j, k;
logical trouve;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* FILL THE TABLE OF POSITIONING POSUIV WHICH ALLOWS TO */
/* PARSE BY COLUMN THE INFERIOR TRIANGULAR PART OF THE */
/* MATRIX IN FORM OF PROFILE */
/* KEYWORDS : */
/* ----------- */
/* RESERVE, MATRIX, PROFILE */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* NISTOC: NUMBER OF COEFFICIENTS IN THE PROFILE */
/* DIMMAT: NUMBER OF LINE OF THE SYMMETRIC SQUARE MATRIX */
/* APOSIT: TABLE OF POSITIONING OF STORAGE TERMS */
/* APOSIT(1,I) CONTAINS THE NUMBER OF TERMES-1 ON LINE */
/* I IN THE PROFILE OF THE MATRIX */
/* APOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF DIAGONAL TERM */
/* OF LINE I */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* POSUIV: POSUIV(K) (WHERE K IS THE INDEX OF STORAGE OF MAT(I,J)) */
/* CONTAINS THE SMALLEST NUMBER IMIN>I OF THE LINE THAT */
/* POSSESSES A TERM MAT(IMIN,J) THAT IS IN THE PROFILE. */
/* IF THERE IS NO TERM MAT(IMIN,J) IN THE PROFILE THEN POSUIV(K)=-1 */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* INITIALIZATIONS */
/* ***********************************************************************
*/
/* Parameter adjustments */
aposit -= 3;
--posuiv;
/* Function Body */
ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
if (ldbg) {
AdvApp2Var_SysBase::mgenmsg_("MMPOSUI", 7L);
}
*iercod = 0;
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
i__1 = *dimmat;
for (i__ = 1; i__ <= i__1; ++i__) {
jmin = i__ - aposit[(i__ << 1) + 1];
i__2 = i__;
for (j = jmin; j <= i__2; ++j) {
imin = i__ + 1;
trouve = FALSE_;
while(! trouve && imin <= *dimmat) {
if (imin - aposit[(imin << 1) + 1] <= j) {
trouve = TRUE_;
} else {
++imin;
}
}
k = aposit[(i__ << 1) + 2] - i__ + j;
if (trouve) {
posuiv[k] = imin;
} else {
posuiv[k] = -1;
}
}
}
goto L9999;
/* ***********************************************************************
*/
/* ERROR PROCESSING */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
L9999:
/* ___ DESALLOCATION, ... */
AdvApp2Var_SysBase::maermsg_("MMPOSUI", iercod, 7L);
if (ldbg) {
AdvApp2Var_SysBase::mgsomsg_("MMPOSUI", 7L);
}
return 0 ;
} /* mmposui_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmresol_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmresol_(integer *hdimen,
integer *gdimen,
integer *hnstoc,
integer *gnstoc,
integer *mnstoc,
doublereal *matsyh,
doublereal *matsyg,
doublereal *vecsyh,
doublereal *vecsyg,
integer *hposit,
integer *hposui,
integer *gposit,
integer *mmposui,
integer *mposit,
doublereal *vecsol,
integer *iercod)
{
integer c__100 = 100;
/* System generated locals */
integer i__1, i__2;
/* Local variables */
logical ldbg;
doublereal* mcho = 0;
integer jmin, jmax, i__, j, k, l;
intptr_t iofv1, iofv2, iofv3, iofv4;
doublereal *v1 = 0, *v2 = 0, *v3 = 0, *v4 = 0;
integer deblig, dimhch;
doublereal* hchole = 0;
intptr_t iofmch, iofmam, iofhch;
doublereal* matsym = 0;
integer ier;
integer aux;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* SOLUTION OF THE SYSTEM */
/* H t(G) V B */
/* = */
/* G 0 L C */
/* KEYWORDS : */
/* ----------- */
/* RESERVE, SOLUTION, SYSTEM, LAGRANGIAN */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* HDIMEN: NOMBER OF LINE (OR COLUMN) OF THE HESSIAN MATRIX */
/* GDIMEN: NOMBER OF LINE OF THE MATRIX OF CONSTRAINTS */
/* HNSTOC: NOMBErS OF TERMS IN THE PROFILE OF HESSIAN MATRIX
*/
/* GNSTOC: NOMBERS OF TERMS IN THE PROFILE OF THE MATRIX OF CONSTRAINTS */
/* MNSTOC: NOMBERS OF TERMS IN THE PROFILE OF THE MATRIX M= G H t(G) */
/* where H IS THE HESSIAN MATRIX AND G IS THE MATRIX OF CONSTRAINTS */
/* MATSYH: TRIANGULAR INFERIOR PART OF THE HESSIAN MATRIX */
/* IN FORM OF PROFILE */
/* MATSYG: MATRIX OF CONSTRAINTS IN FORM OF PROFILE */
/* VECSYH: VECTOR OF THE SECOND MEMBER ASSOCIATED TO MATSYH */
/* VECSYG: VECTOR OF THE SECOND MEMBER ASSOCIATED TO MATSYG */
/* HPOSIT: TABLE OF POSITIONING OF THE HESSIAN MATRIX */
/* HPOSIT(1,I) CONTAINS THE NUMBER OF TERMS -1 */
/* WHICH ARE IN THE PROFILE AT LINE I */
/* HPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF TERM */
/* DIAGONAL OF THE MATRIX AT LINE I */
/* HPOSUI: TABLE ALLOWING TO PARSE THE HESSIAN MATRIX BY COLUMN */
/* IN FORM OF PROFILE */
/* HPOSUI(K) CONTAINS THE NUMBER OF LINE IMIN FOLLOWING THE CURRENT LINE*/
/* I WHERE H(I,J)=MATSYH(K) AS IT EXISTS IN THE */
/* SAME COLUMN J A TERM IN THE PROFILE OF LINE IMIN */
/* IF SUCH TERM DOES NOT EXIST IMIN=-1 */
/* GPOSIT: TABLE OF POSITIONING OF THE MATRIX OF CONSTRAINTS */
/* GPOSIT(1,I) CONTAINS THE NUMBER OF TERMS OF LINE I */
/* WHICH ARE IN THE PROFILE */
/* GPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF THE LAST TERM */
/* OF LINE I WHICH IS IN THE PROFILE */
/* GPOSIT(3,I) CONTAINS THE NUMBER OF COLUMN CORRESPONDING */
/* TO THE FIRST TERM OF LINE I WHICH IS IN THE PROFILE */
/* MMPOSUI, MPOSIT: SAME STRUCTURE AS HPOSUI, BUT FOR MATRIX */
/* M=G H t(G) */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* VECSOL: VECTOR SOLUTION V OF THE SYSTEM */
/* IERCOD: ERROR CODE */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* INITIALISATIONS */
/* ***********************************************************************
*/
/* Parameter adjustments */
--vecsol;
hposit -= 3;
--vecsyh;
--hposui;
--matsyh;
--matsyg;
--vecsyg;
gposit -= 4;
--mmposui;
mposit -= 3;
/* Function Body */
ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
if (ldbg) {
AdvApp2Var_SysBase::mgenmsg_("MMRESOL", 7L);
}
*iercod = 0;
iofhch = 0;
iofv1 = 0;
iofv2 = 0;
iofv3 = 0;
iofv4 = 0;
iofmam = 0;
iofmch = 0;
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
/* Dynamic allocation */
AdvApp2Var_SysBase anAdvApp2Var_SysBase;
anAdvApp2Var_SysBase.macrar8_(hdimen, &c__100, v1, &iofv1, &ier);
if (ier > 0) {
goto L9102;
}
dimhch = hposit[(*hdimen << 1) + 2];
anAdvApp2Var_SysBase.macrar8_(&dimhch, &c__100, hchole, &iofhch, &ier);
if (ier > 0) {
goto L9102;
}
/* solution of system 1 H V1 = b */
/* where H=MATSYH and b=VECSYH */
mmchole_(hnstoc, hdimen, &matsyh[1], &hposit[3], &hposui[1], &hchole[
iofhch], &ier);
if (ier > 0) {
goto L9101;
}
mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1], &vecsyh[
1], &v1[iofv1], &ier);
if (ier > 0) {
goto L9102;
}
/* Case when there are constraints */
if (*gdimen > 0) {
/* Calculate the vector of the second member V2=G H(-1) b -c = G v1-c */
/* of system of unknown Lagrangian vector MULTIP */
/* where G=MATSYG */
/* c=VECSYG */
anAdvApp2Var_SysBase.macrar8_(gdimen, &c__100, v2, &iofv2, &ier);
if (ier > 0) {
goto L9102;
}
anAdvApp2Var_SysBase.macrar8_(hdimen, &c__100, v3, &iofv3, &ier);
if (ier > 0) {
goto L9102;
}
anAdvApp2Var_SysBase.macrar8_(gdimen, &c__100, v4, &iofv4, &ier);
if (ier > 0) {
goto L9102;
}
anAdvApp2Var_SysBase.macrar8_(mnstoc, &c__100, matsym, &iofmam, &ier);
if (ier > 0) {
goto L9102;
}
deblig = 1;
mmatvec_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v1[iofv1], &
deblig, &v2[iofv2], &ier);
if (ier > 0) {
goto L9101;
}
i__1 = *gdimen;
for (i__ = 1; i__ <= i__1; ++i__) {
v2[i__ + iofv2 - 1] -= vecsyg[i__];
}
/* Calculate the matrix M= G H(-1) t(G) */
/* RESOL DU SYST 2 : H qi = gi */
/* where is a vector column of t(G) */
/* qi=v3 */
/* then calculate G qi */
/* then construct M in form of profile */
i__1 = *gdimen;
for (i__ = 1; i__ <= i__1; ++i__) {
AdvApp2Var_SysBase::mvriraz_(hdimen, &v1[iofv1]);
AdvApp2Var_SysBase::mvriraz_(hdimen, &v3[iofv3]);
AdvApp2Var_SysBase::mvriraz_(gdimen, &v4[iofv4]);
jmin = gposit[i__ * 3 + 3];
jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
i__2 = jmax;
for (j = jmin; j <= i__2; ++j) {
k = j + aux;
v1[j + iofv1 - 1] = matsyg[k];
}
mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1],
&v1[iofv1], &v3[iofv3], &ier);
if (ier > 0) {
goto L9101;
}
deblig = i__;
mmatvec_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v3[
iofv3], &deblig, &v4[iofv4], &ier);
if (ier > 0) {
goto L9101;
}
k = mposit[(i__ << 1) + 2];
matsym[k + iofmam - 1] = v4[i__ + iofv4 - 1];
while(mmposui[k] > 0) {
l = mmposui[k];
k = mposit[(l << 1) + 2] - l + i__;
matsym[k + iofmam - 1] = v4[l + iofv4 - 1];
}
}
/* SOLVE SYST 3 M L = V2 */
/* WITH L=V4 */
AdvApp2Var_SysBase::mvriraz_(gdimen, &v4[iofv4]);
anAdvApp2Var_SysBase.macrar8_(mnstoc, &c__100, mcho, &iofmch, &ier);
if (ier > 0) {
goto L9102;
}
mmchole_(mnstoc, gdimen, &matsym[iofmam], &mposit[3], &mmposui[1], &
mcho[iofmch], &ier);
if (ier > 0) {
goto L9101;
}
mmrslss_(mnstoc, gdimen, &mcho[iofmch], &mposit[3], &mmposui[1], &v2[
iofv2], &v4[iofv4], &ier);
if (ier > 0) {
goto L9102;
}
/* CALCULATE THE VECTOR OF THE SECOND MEMBER OF THE SYSTEM Hx = b - t(G) L
*/
/* = V1 */
AdvApp2Var_SysBase::mvriraz_(hdimen, &v1[iofv1]);
mmtmave_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v4[iofv4], &
v1[iofv1], &ier);
if (ier > 0) {
goto L9101;
}
i__1 = *hdimen;
for (i__ = 1; i__ <= i__1; ++i__) {
v1[i__ + iofv1 - 1] = vecsyh[i__] - v1[i__ + iofv1 - 1];
}
/* RESOL SYST 4 Hx = b - t(G) L */
mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1], &v1[
iofv1], &vecsol[1], &ier);
if (ier > 0) {
goto L9102;
}
} else {
i__1 = *hdimen;
for (i__ = 1; i__ <= i__1; ++i__) {
vecsol[i__] = v1[i__ + iofv1 - 1];
}
}
goto L9999;
/* ***********************************************************************
*/
/* PROCESSING OF ERRORS */
/* ***********************************************************************
*/
L9101:
*iercod = 1;
goto L9999;
L9102:
AdvApp2Var_SysBase::mswrdbg_("MMRESOL : PROBLEM WITH DIMMAT", 30L);
*iercod = 2;
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
L9999:
/* ___ DESALLOCATION, ... */
anAdvApp2Var_SysBase.macrdr8_(hdimen, &c__100, v1, &iofv1, &ier);
if (*iercod == 0 && ier > 0) {
*iercod = 3;
}
anAdvApp2Var_SysBase.macrdr8_(&dimhch, &c__100, hchole, &iofhch, &ier);
if (*iercod == 0 && ier > 0) {
*iercod = 3;
}
anAdvApp2Var_SysBase.macrdr8_(gdimen, &c__100, v2, &iofv2, &ier);
if (*iercod == 0 && ier > 0) {
*iercod = 3;
}
anAdvApp2Var_SysBase.macrdr8_(hdimen, &c__100, v3, &iofv3, &ier);
if (*iercod == 0 && ier > 0) {
*iercod = 3;
}
anAdvApp2Var_SysBase.macrdr8_(gdimen, &c__100, v4, &iofv4, &ier);
if (*iercod == 0 && ier > 0) {
*iercod = 3;
}
anAdvApp2Var_SysBase.macrdr8_(mnstoc, &c__100, matsym, &iofmam, &ier);
if (*iercod == 0 && ier > 0) {
*iercod = 3;
}
anAdvApp2Var_SysBase.macrdr8_(mnstoc, &c__100, mcho, &iofmch, &ier);
if (*iercod == 0 && ier > 0) {
*iercod = 3;
}
AdvApp2Var_SysBase::maermsg_("MMRESOL", iercod, 7L);
if (ldbg) {
AdvApp2Var_SysBase::mgsomsg_("MMRESOL", 7L);
}
return 0 ;
} /* mmresol_ */
//=======================================================================
//function : mmrslss_
//purpose :
//=======================================================================
int mmrslss_(integer *,//mxcoef,
integer *dimens,
doublereal *smatri,
integer *sposit,
integer *posuiv,
doublereal *mscnmbr,
doublereal *soluti,
integer *iercod)
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
logical ldbg;
integer i__, j;
doublereal somme;
integer pointe, ptcour;
/* ***********************************************************************
*/
/* FuNCTION : */
/* ---------- T */
/* Solves linear system SS x = b where S is a */
/* triangular lower matrix given in form of profile */
/* KEYWORDS : */
/* ----------- */
/* RESERVE, MATRICE_PROFILE, RESOLUTION, CHOLESKI */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* MXCOEF : Maximum number of non-null coefficient in the matrix */
/* DIMENS : Dimension of the matrix */
/* SMATRI(MXCOEF) : Values of coefficients of the matrix */
/* SPOSIT(2,DIMENS): */
/* SPOSIT(1,*) : Distance diagonal-extremity of the line */
/* SPOSIT(2,*) : Position of diagonal terms in AMATRI */
/* POSUIV(MXCOEF): first line inferior not out of profile */
/* MSCNMBR(DIMENS): Vector second member of the equation */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* SOLUTI(NDIMEN) : Result vector */
/* IERCOD : Error code 0 : ok */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* T */
/* SS is the decomposition of choleski of a symmetric matrix */
/* defined postive, that can result from routine MMCHOLE. */
/* For a full matrix it is possible to use MRSLMSC */
/* LEVEL OF DEBUG = 4 */
/* > */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* INITIALISATIONS */
/* ***********************************************************************
*/
/* Parameter adjustments */
--posuiv;
--smatri;
--soluti;
--mscnmbr;
sposit -= 3;
/* Function Body */
ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
if (ldbg) {
AdvApp2Var_SysBase::mgenmsg_("MMRSLSS", 7L);
}
*iercod = 0;
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
/* ----- Solution of Sw = b */
i__1 = *dimens;
for (i__ = 1; i__ <= i__1; ++i__) {
pointe = sposit[(i__ << 1) + 2];
somme = 0.;
i__2 = i__ - 1;
for (j = i__ - sposit[(i__ << 1) + 1]; j <= i__2; ++j) {
somme += smatri[pointe - (i__ - j)] * soluti[j];
}
soluti[i__] = (mscnmbr[i__] - somme) / smatri[pointe];
}
/* T */
/* ----- Solution of S u = w */
for (i__ = *dimens; i__ >= 1; --i__) {
pointe = sposit[(i__ << 1) + 2];
j = posuiv[pointe];
somme = 0.;
while(j > 0) {
ptcour = sposit[(j << 1) + 2] - (j - i__);
somme += smatri[ptcour] * soluti[j];
j = posuiv[ptcour];
}
soluti[i__] = (soluti[i__] - somme) / smatri[pointe];
}
goto L9999;
/* ***********************************************************************
*/
/* ERROR PROCESSING */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* RETURN PROGRAM CALLING */
/* ***********************************************************************
*/
L9999:
AdvApp2Var_SysBase::maermsg_("MMRSLSS", iercod, 7L);
if (ldbg) {
AdvApp2Var_SysBase::mgsomsg_("MMRSLSS", 7L);
}
return 0 ;
} /* mmrslss_ */
//=======================================================================
//function : mmrslw_
//purpose :
//=======================================================================
int mmrslw_(integer *normax,
integer *nordre,
integer *ndimen,
doublereal *epspiv,
doublereal *abmatr,
doublereal *xmatri,
integer *iercod)
{
/* System generated locals */
integer abmatr_dim1, abmatr_offset, xmatri_dim1, xmatri_offset, i__1,
i__2, i__3;
doublereal d__1;
/* Local variables */
integer kpiv;
doublereal pivot;
integer ii, jj, kk;
doublereal akj;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Solution of a linear system A.x = B of N equations to N */
/* unknown by Gauss method (partial pivot) or : */
/* A is matrix NORDRE * NORDRE, */
/* B is matrix NORDRE (lines) * NDIMEN (columns), */
/* x is matrix NORDRE (lines) * NDIMEN (columns). */
/* In this program, A and B are stored in matrix ABMATR */
/* the lines and columns which of were inverted. ABMATR(k,j) is */
/* term A(j,k) if k <= NORDRE, B(j,k-NORDRE) otherwise (see example). */
/* KEYWORDS : */
/* ----------- */
/* TOUS, MATH_ACCES::EQUATION&, MATRICE&, RESOLUTION, GAUSS, &SOLUTION */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NORMAX : Max size of the first index of XMATRI. This argument */
/* serves only for the declaration of dimension of XMATRI and should be */
/* above or equal to NORDRE. */
/* NORDRE : Order of the matrix i.e. number of equations and */
/* unknown quantities of the linear system to be solved. */
/* NDIMEN : Number of the second member. */
/* EPSPIV : Minimal value of a pivot. If during the calculation */
/* the absolute value of the pivot is below EPSPIV, the */
/* system of equations is declared singular. EPSPIV should */
/* be a "small" real. */
/* ABMATR(NORDRE+NDIMEN,NORDRE) : Auxiliary matrix containing */
/* matrix A and matrix B. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* XMATRI : Matrix containing NORDRE*NDIMEN solutions. */
/* IERCOD=0 shows that all solutions are calculated. */
/* IERCOD=1 shows that the matrix is of lower rank than NORDRE */
/* (the system is singular). */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ATTENTION : the indices of line and column are inverted */
/* compared to usual indices. */
/* System : */
/* a1*x + b1*y = c1 */
/* a2*x + b2*y = c2 */
/* should be represented by matrix ABMATR : */
/* ABMATR(1,1) = a1 ABMATR(1,2) = a2 */
/* ABMATR(2,1) = b1 ABMATR(2,2) = b2 */
/* ABMATR(3,1) = c1 ABMATR(3,2) = c2 */
/* To solve this system, it is necessary to set : */
/* NORDRE = 2 (there are 2 equations with 2 unknown values), */
/* NDIMEN = 1 (there is only one second member), */
/* any NORMAX can be taken >= NORDRE. */
/* To use this routine, it is recommended to use one of */
/* interfaces : MMRSLWI or MMMRSLWD. */
/* > */
/* **********************************************************************
*/
/* Name of the routine */
/* INTEGER IBB,MNFNDEB */
/* IBB=MNFNDEB() */
/* IF (IBB.GE.2) CALL MGENMSG(NOMPR) */
/* Parameter adjustments */
xmatri_dim1 = *normax;
xmatri_offset = xmatri_dim1 + 1;
xmatri -= xmatri_offset;
abmatr_dim1 = *nordre + *ndimen;
abmatr_offset = abmatr_dim1 + 1;
abmatr -= abmatr_offset;
/* Function Body */
*iercod = 0;
/* *********************************************************************
*/
/* Triangulation of matrix ABMATR. */
/* *********************************************************************
*/
i__1 = *nordre;
for (kk = 1; kk <= i__1; ++kk) {
/* ---------- Find max pivot in column KK. ------------
--- */
pivot = *epspiv;
kpiv = 0;
i__2 = *nordre;
for (jj = kk; jj <= i__2; ++jj) {
akj = (d__1 = abmatr[kk + jj * abmatr_dim1], advapp_abs(d__1));
if (akj > pivot) {
pivot = akj;
kpiv = jj;
}
/* L100: */
}
if (kpiv == 0) {
goto L9900;
}
/* --------- Swapping of line KPIV with line KK. ------
--- */
if (kpiv != kk) {
i__2 = *nordre + *ndimen;
for (jj = kk; jj <= i__2; ++jj) {
akj = abmatr[jj + kk * abmatr_dim1];
abmatr[jj + kk * abmatr_dim1] = abmatr[jj + kpiv *
abmatr_dim1];
abmatr[jj + kpiv * abmatr_dim1] = akj;
/* L200: */
}
}
/* ---------- Removal and triangularization. -----------
--- */
pivot = -abmatr[kk + kk * abmatr_dim1];
i__2 = *nordre;
for (ii = kk + 1; ii <= i__2; ++ii) {
akj = abmatr[kk + ii * abmatr_dim1] / pivot;
i__3 = *nordre + *ndimen;
for (jj = kk + 1; jj <= i__3; ++jj) {
abmatr[jj + ii * abmatr_dim1] += akj * abmatr[jj + kk *
abmatr_dim1];
/* L400: */
}
/* L300: */
}
/* L1000: */
}
/* *********************************************************************
*/
/* Solution of the system of triangular equations. */
/* Matrix ABMATR(NORDRE+JJ,II), contains second members */
/* of the system for 1<=j<=NDIMEN and 1<=i<=NORDRE. */
/* *********************************************************************
*/
/* ---------------- Calculation of solutions by ascending. -----------------
*/
for (kk = *nordre; kk >= 1; --kk) {
pivot = abmatr[kk + kk * abmatr_dim1];
i__1 = *ndimen;
for (ii = 1; ii <= i__1; ++ii) {
akj = abmatr[ii + *nordre + kk * abmatr_dim1];
i__2 = *nordre;
for (jj = kk + 1; jj <= i__2; ++jj) {
akj -= abmatr[jj + kk * abmatr_dim1] * xmatri[jj + ii *
xmatri_dim1];
/* L800: */
}
xmatri[kk + ii * xmatri_dim1] = akj / pivot;
/* L700: */
}
/* L600: */
}
goto L9999;
/* ------If the absolute value of a pivot is smaller than -------- */
/* ---------- EPSPIV: return the code of error. ------------
*/
L9900:
*iercod = 1;
L9999:
if (*iercod > 0) {
AdvApp2Var_SysBase::maermsg_("MMRSLW ", iercod, 7L);
}
/* IF (IBB.GE.2) CALL MGSOMSG(NOMPR) */
return 0 ;
} /* mmrslw_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmmrslwd_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmmrslwd_(integer *normax,
integer *nordre,
integer *ndim,
doublereal *amat,
doublereal *bmat,
doublereal *epspiv,
doublereal *aaux,
doublereal *xmat,
integer *iercod)
{
/* System generated locals */
integer amat_dim1, amat_offset, bmat_dim1, bmat_offset, xmat_dim1,
xmat_offset, aaux_dim1, aaux_offset, i__1, i__2;
/* Local variables */
integer i__, j;
integer ibb;
/* IMPLICIT DOUBLE PRECISION (A-H,O-Z) */
/* IMPLICIT INTEGER (I-N) */
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Solution of a linear system by Gauss method where */
/* the second member is a table of vectors. Method of partial pivot. */
/* KEYWORDS : */
/* ----------- */
/* ALL, MATH_ACCES :: */
/* SYSTEME&,EQUATION&, RESOLUTION,GAUSS ,&VECTEUR */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NORMAX : Max. Dimension of AMAT. */
/* NORDRE : Order of the matrix. */
/* NDIM : Number of columns of BMAT and XMAT. */
/* AMAT(NORMAX,NORDRE) : The processed matrix. */
/* BMAT(NORMAX,NDIM) : The matrix of second member. */
/* XMAT(NORMAX,NDIM) : The matrix of solutions. */
/* EPSPIV : Min value of a pivot. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* AAUX(NORDRE+NDIM,NORDRE) : Auxiliary matrix. */
/* XMAT(NORMAX,NDIM) : Matrix of solutions. */
/* IERCOD=0 shows that solutions in XMAT are valid. */
/* IERCOD=1 shows that matrix AMAT is of lower rank than NORDRE. */
/* COMMONS USED : */
/* ---------------- */
/* .Neant. */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Type Name */
/* MAERMSG MGENMSG MGSOMSG */
/* MMRSLW I*4 MNFNDEB */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ATTENTION : lines and columns are located in usual order : */
/* 1st index = index line */
/* 2nd index = index column */
/* Example, the system : */
/* a1*x + b1*y = c1 */
/* a2*x + b2*y = c2 */
/* is represented by matrix AMAT : */
/* AMAT(1,1) = a1 AMAT(2,1) = a2 */
/* AMAT(1,2) = b1 AMAT(2,2) = b2 */
/* The first index is the index of line, the second index */
/* is the index of columns (Compare with MMRSLWI which is faster). */
/* > */
/* **********************************************************************
*/
/* Name of the routine */
/* Parameter adjustments */
amat_dim1 = *normax;
amat_offset = amat_dim1 + 1;
amat -= amat_offset;
xmat_dim1 = *normax;
xmat_offset = xmat_dim1 + 1;
xmat -= xmat_offset;
aaux_dim1 = *nordre + *ndim;
aaux_offset = aaux_dim1 + 1;
aaux -= aaux_offset;
bmat_dim1 = *normax;
bmat_offset = bmat_dim1 + 1;
bmat -= bmat_offset;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 3) {
AdvApp2Var_SysBase::mgenmsg_("MMMRSLW", 7L);
}
/* Initialization of the auxiliary matrix. */
i__1 = *nordre;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = *nordre;
for (j = 1; j <= i__2; ++j) {
aaux[j + i__ * aaux_dim1] = amat[i__ + j * amat_dim1];
/* L200: */
}
/* L100: */
}
/* Second member. */
i__1 = *nordre;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = *ndim;
for (j = 1; j <= i__2; ++j) {
aaux[j + *nordre + i__ * aaux_dim1] = bmat[i__ + j * bmat_dim1];
/* L400: */
}
/* L300: */
}
/* Solution of the system of equations. */
mmrslw_(normax, nordre, ndim, epspiv, &aaux[aaux_offset], &xmat[
xmat_offset], iercod);
if (*iercod != 0) {
AdvApp2Var_SysBase::maermsg_("MMMRSLW", iercod, 7L);
}
if (ibb >= 3) {
AdvApp2Var_SysBase::mgsomsg_("MMMRSLW", 7L);
}
return 0 ;
} /* mmmrslwd_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmrtptt_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmrtptt_(integer *ndglgd,
doublereal *rtlegd)
{
integer ideb, nmod2, nsur2, ilong, ibb;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Extracts from Common LDGRTL the STRICTLY positive roots of the */
/* Legendre polynom of degree NDGLGD, for 2 <= NDGLGD <= 61. */
/* KEYWORDS : */
/* ----------- */
/* TOUS, AB_SPECIFI::COMMON&, EXTRACTION, &RACINE, &LEGENDRE. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDGLGD : Mathematic degree of Legendre polynom. */
/* This degree should be above or equal to 2 and */
/* below or equal to 61. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* RTLEGD : The table of strictly positive roots of */
/* Legendre polynom of degree NDGLGD. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* ATTENTION: the condition on NDEGRE ( 2 <= NDEGRE <= 61) is not */
/* tested. The caller should make the test. */
/* > */
/* **********************************************************************
*/
/* Nome of the routine */
/* Common MLGDRTL: */
/* This common includes POSITIVE roots of Legendre polynoms */
/* AND the weight of Gauss quadrature formulas on all */
/* POSITIVE roots of Legendre polynoms. */
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* The common of Legendre roots. */
/* KEYWORDS : */
/* ----------- */
/* BASE LEGENDRE */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* ROOTAB : Table of all rotts of Legendre polynoms */
/* between [0,1]. They are ranked for degrees increasing from 2 to 61. */
/* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
/* The address is the same. */
/* HI0TAB : Table of Legendre interpolators for root x=0 */
/* the polynoms of UNEVEN degree. */
/* RTLTB0 : Table of Li(uk) where uk are roots of a */
/* Legendre polynom of EVEN degree. */
/* RTLTB1 : Table of Li(uk) where uk are roots of a */
/* Legendre polynom of UNEVEN degree. */
/************************************************************************
*****/
/* Parameter adjustments */
--rtlegd;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 3) {
AdvApp2Var_SysBase::mgenmsg_("MMRTPTT", 7L);
}
if (*ndglgd < 2) {
goto L9999;
}
nsur2 = *ndglgd / 2;
nmod2 = *ndglgd % 2;
ilong = nsur2 << 3;
ideb = nsur2 * (nsur2 - 1) / 2 + 1;
AdvApp2Var_SysBase::mcrfill_(&ilong,
&mlgdrtl_.rootab[ideb + nmod2 * 465 - 1],
&rtlegd[1]);
/* ----------------------------- The end --------------------------------
*/
L9999:
if (ibb >= 3) {
AdvApp2Var_SysBase::mgsomsg_("MMRTPTT", 7L);
}
return 0;
} /* mmrtptt_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmsrre2_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmsrre2_(doublereal *tparam,
integer *nbrval,
doublereal *tablev,
doublereal *epsil,
integer *numint,
integer *itypen,
integer *iercod)
{
/* System generated locals */
doublereal d__1;
/* Local variables */
integer ideb, ifin, imil, ibb;
/* ***********************************************************************
*/
/* FUNCTION : */
/* -------- */
/* Find the interval corresponding to a valueb given in */
/* increasing order of real numbers with double precision. */
/* KEYWORDS : */
/* --------- */
/* TOUS,MATH_ACCES::TABLEAU&,POINT&,CORRESPONDANCE,&RANG */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* TPARAM : Value to be tested. */
/* NBRVAL : Size of TABLEV */
/* TABLEV : Table of reals. */
/* EPSIL : Epsilon of precision */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* NUMINT : Number of the interval (between 1 and NBRVAL-1). */
/* ITYPEN : = 0 TPARAM is inside the interval NUMINT */
/* = 1 : TPARAM corresponds to the lower limit of */
/* the provided interval. */
/* = 2 : TPARAM corresponds to the upper limit of */
/* the provided interval. */
/* IERCOD : Error code. */
/* = 0 : OK */
/* = 1 : TABLEV does not contain enough elements. */
/* = 2 : TPARAM out of limits of TABLEV. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* --------------------------------- */
/* There are NBRVAL values in TABLEV which stands for NBRVAL-1 intervals. */
/* One searches the interval containing TPARAM by */
/* dichotomy. Complexity of the algorithm : Log(n)/Log(2).(RBD). */
/* > */
/* ***********************************************************************
*/
/* Initialisations */
/* Parameter adjustments */
--tablev;
/* Function Body */
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 6) {
AdvApp2Var_SysBase::mgenmsg_("MMSRRE2", 7L);
}
*iercod = 0;
*numint = 0;
*itypen = 0;
ideb = 1;
ifin = *nbrval;
/* TABLEV should contain at least two values */
if (*nbrval < 2) {
*iercod = 1;
goto L9999;
}
/* TPARAM should be between extreme limits of TABLEV. */
if (*tparam < tablev[1] || *tparam > tablev[*nbrval]) {
*iercod = 2;
goto L9999;
}
/* ----------------------- SEARCH OF THE INTERVAL --------------------
*/
L1000:
/* Test end of loop (found). */
if (ideb + 1 == ifin) {
*numint = ideb;
goto L2000;
}
/* Find by dichotomy on increasing values of TABLEV. */
imil = (ideb + ifin) / 2;
if (*tparam >= tablev[ideb] && *tparam <= tablev[imil]) {
ifin = imil;
} else {
ideb = imil;
}
goto L1000;
/* -------------- TEST IF TPARAM IS NOT A VALUE --------- */
/* ------------------------OF TABLEV UP TO EPSIL ----------------------
*/
L2000:
if ((d__1 = *tparam - tablev[ideb], advapp_abs(d__1)) < *epsil) {
*itypen = 1;
goto L9999;
}
if ((d__1 = *tparam - tablev[ifin], advapp_abs(d__1)) < *epsil) {
*itypen = 2;
goto L9999;
}
/* --------------------------- THE END ----------------------------------
*/
L9999:
if (*iercod > 0) {
AdvApp2Var_SysBase::maermsg_("MMSRRE2", iercod, 7L);
}
if (ibb >= 6) {
AdvApp2Var_SysBase::mgsomsg_("MMSRRE2", 7L);
}
return 0 ;
} /* mmsrre2_ */
//=======================================================================
//function : mmtmave_
//purpose :
//=======================================================================
int mmtmave_(integer *nligne,
integer *ncolon,
integer *gposit,
integer *,//gnstoc,
doublereal *gmatri,
doublereal *vecin,
doublereal *vecout,
integer *iercod)
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
logical ldbg;
integer imin, imax, i__, j, k;
doublereal somme;
integer aux;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* t */
/* CREATES PRODUCT G V */
/* WHERE THE MATRIX IS IN FORM OF PROFILE */
/* KEYWORDS : */
/* ----------- */
/* RESERVE, PRODUCT, MATRIX, PROFILE, VECTOR */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* NLIGNE : NUMBER OF LINE OF THE MATRIX */
/* NCOLON : NOMBER OF COLUMN OF THE MATRIX */
/* GPOSIT: TABLE OF POSITIONING OF TERMS OF STORAGE */
/* GPOSIT(1,I) CONTAINS THE NUMBER of TERMS-1 ON LINE */
/* I IN THE PROFILE OF THE MATRIX */
/* GPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF THE DIAGONAL TERM*/
/* OF LINE I */
/* GPOSIT(3,I) CONTAINS THE INDEX COLUMN OF THE FIRST TERM OF */
/* PROFILE OF LINE I */
/* GNSTOC : NOMBER OF TERM IN THE PROFILE OF GMATRI */
/* GMATRI : MATRIX OF CONSTRAINTS IN FORM OF PROFILE */
/* VECIN : INPUT VECTOR */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* VECOUT : VECTOR PRODUCT */
/* IERCOD : ERROR CODE */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* INITIALISATIONS */
/* ***********************************************************************
*/
/* Parameter adjustments */
--vecin;
gposit -= 4;
--vecout;
--gmatri;
/* Function Body */
ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
if (ldbg) {
AdvApp2Var_SysBase::mgenmsg_("MMTMAVE", 7L);
}
*iercod = 0;
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
i__1 = *ncolon;
for (i__ = 1; i__ <= i__1; ++i__) {
somme = 0.;
i__2 = *nligne;
for (j = 1; j <= i__2; ++j) {
imin = gposit[j * 3 + 3];
imax = gposit[j * 3 + 1] + gposit[j * 3 + 3] - 1;
aux = gposit[j * 3 + 2] - gposit[j * 3 + 1] - imin + 1;
if (imin <= i__ && i__ <= imax) {
k = i__ + aux;
somme += gmatri[k] * vecin[j];
}
}
vecout[i__] = somme;
}
goto L9999;
/* ***********************************************************************
*/
/* ERROR PROCESSING */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
L9999:
/* ___ DESALLOCATION, ... */
AdvApp2Var_SysBase::maermsg_("MMTMAVE", iercod, 7L);
if (ldbg) {
AdvApp2Var_SysBase::mgsomsg_("MMTMAVE", 7L);
}
return 0 ;
} /* mmtmave_ */
//=======================================================================
//function : mmtrpj0_
//purpose :
//=======================================================================
int mmtrpj0_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *epsi3d,
doublereal *crvlgd,
doublereal *ycvmax,
doublereal *epstrc,
integer *ncfnew)
{
/* System generated locals */
integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer ncut, i__;
doublereal bidon, error;
integer nd;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Lowers the degree of a curve defined on (-1,1) in the direction of */
/* Legendre with a given precision. */
/* KEYWORDS : */
/* ----------- */
/* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX : Max Nb of coeff. of the curve (dimensioning). */
/* NDIMEN : Dimension of the space. */
/* NCOEFF : Degree +1 of the polynom. */
/* EPSI3D : Precision required for the approximation. */
/* CRVLGD : The curve the degree which of it is required to lower. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* EPSTRC : Precision of the approximation. */
/* NCFNEW : Degree +1 of the resulting polynom. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* ------- Minimum degree that can be attained : Stop at 1 (RBD) ---------
*/
/* Parameter adjustments */
--ycvmax;
crvlgd_dim1 = *ncofmx;
crvlgd_offset = crvlgd_dim1 + 1;
crvlgd -= crvlgd_offset;
/* Function Body */
*ncfnew = 1;
/* ------------------- Init for error calculation -----------------------
*/
i__1 = *ndimen;
for (i__ = 1; i__ <= i__1; ++i__) {
ycvmax[i__] = 0.;
/* L100: */
}
*epstrc = 0.;
error = 0.;
/* Cutting of coefficients. */
ncut = 2;
/* ------ Loop on the series of Legendre :NCOEFF --> 2 (RBD) -----------
*/
i__1 = ncut;
for (i__ = *ncoeff; i__ >= i__1; --i__) {
/* Factor of renormalization. */
bidon = ((i__ - 1) * 2. + 1.) / 2.;
bidon = sqrt(bidon);
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
bidon;
/* L310: */
}
/* Cutting is stopped if the norm becomes too great. */
error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
if (error > *epsi3d) {
*ncfnew = i__;
goto L9999;
}
/* --- Max error cumulee when the I-th coeff is removed. */
*epstrc = error;
/* L300: */
}
/* --------------------------------- End --------------------------------
*/
L9999:
return 0;
} /* mmtrpj0_ */
//=======================================================================
//function : mmtrpj2_
//purpose :
//=======================================================================
int mmtrpj2_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *epsi3d,
doublereal *crvlgd,
doublereal *ycvmax,
doublereal *epstrc,
integer *ncfnew)
{
/* Initialized data */
static doublereal xmaxj[57] = { .9682458365518542212948163499456,
.986013297183269340427888048593603,
1.07810420343739860362585159028115,
1.17325804490920057010925920756025,
1.26476561266905634732910520370741,
1.35169950227289626684434056681946,
1.43424378958284137759129885012494,
1.51281316274895465689402798226634,
1.5878364329591908800533936587012,
1.65970112228228167018443636171226,
1.72874345388622461848433443013543,
1.7952515611463877544077632304216,
1.85947199025328260370244491818047,
1.92161634324190018916351663207101,
1.98186713586472025397859895825157,
2.04038269834980146276967984252188,
2.09730119173852573441223706382076,
2.15274387655763462685970799663412,
2.20681777186342079455059961912859,
2.25961782459354604684402726624239,
2.31122868752403808176824020121524,
2.36172618435386566570998793688131,
2.41117852396114589446497298177554,
2.45964731268663657873849811095449,
2.50718840313973523778244737914028,
2.55385260994795361951813645784034,
2.59968631659221867834697883938297,
2.64473199258285846332860663371298,
2.68902863641518586789566216064557,
2.73261215675199397407027673053895,
2.77551570192374483822124304745691,
2.8177699459714315371037628127545,
2.85940333797200948896046563785957,
2.90044232019793636101516293333324,
2.94091151970640874812265419871976,
2.98083391718088702956696303389061,
3.02023099621926980436221568258656,
3.05912287574998661724731962377847,
3.09752842783622025614245706196447,
3.13546538278134559341444834866301,
3.17295042316122606504398054547289,
3.2099992681699613513775259670214,
3.24662674946606137764916854570219,
3.28284687953866689817670991319787,
3.31867291347259485044591136879087,
3.35411740487202127264475726990106,
3.38919225660177218727305224515862,
3.42390876691942143189170489271753,
3.45827767149820230182596660024454,
3.49230918177808483937957161007792,
3.5260130200285724149540352829756,
3.55939845146044235497103883695448,
3.59247431368364585025958062194665,
3.62524904377393592090180712976368,
3.65773070318071087226169680450936,
3.68992700068237648299565823810245,
3.72184531357268220291630708234186 };
/* System generated locals */
integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer ncut, i__;
doublereal bidon, error;
integer ia, nd;
doublereal bid, eps1;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Lower the degree of a curve defined on (-1,1) in the direction of */
/* Legendre with a given precision. */
/* KEYWORDS : */
/* ----------- */
/* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
/* NDIMEN : Dimension of the space. */
/* NCOEFF : Degree +1 of the polynom. */
/* EPSI3D : Precision required for the approximation. */
/* CRVLGD : The curve the degree which of will be lowered. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* YCVMAX : Auxiliary table (error max on each dimension).
*/
/* EPSTRC : Precision of the approximation. */
/* NCFNEW : Degree +1 of the resulting polynom. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* Parameter adjustments */
--ycvmax;
crvlgd_dim1 = *ncofmx;
crvlgd_offset = crvlgd_dim1 + 1;
crvlgd -= crvlgd_offset;
/* Function Body */
/* Minimum degree that can be reached : Stop at IA (RBD). -------------
*/
ia = 2;
*ncfnew = ia;
/* Init for calculation of error. */
i__1 = *ndimen;
for (i__ = 1; i__ <= i__1; ++i__) {
ycvmax[i__] = 0.;
/* L100: */
}
*epstrc = 0.;
error = 0.;
/* Cutting of coefficients. */
ncut = ia + 1;
/* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
*/
i__1 = ncut;
for (i__ = *ncoeff; i__ >= i__1; --i__) {
/* Factor of renormalization. */
bidon = xmaxj[i__ - ncut];
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
bidon;
/* L310: */
}
/* One stops to cut if the norm becomes too great. */
error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
if (error > *epsi3d) {
*ncfnew = i__;
goto L400;
}
/* --- Max error cumulated when the I-th coeff is removed. */
*epstrc = error;
/* L300: */
}
/* ------- Cutting of zero coeffs of interpolation (RBD) -------
*/
L400:
if (*ncfnew == ia) {
AdvApp2Var_MathBase::mmeps1_(&eps1);
for (i__ = ia; i__ >= 2; --i__) {
bid = 0.;
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
/* L600: */
}
if (bid > eps1) {
*ncfnew = i__;
goto L9999;
}
/* L500: */
}
/* --- If all coeffs can be removed, this is a point. */
*ncfnew = 1;
}
/* --------------------------------- End --------------------------------
*/
L9999:
return 0;
} /* mmtrpj2_ */
//=======================================================================
//function : mmtrpj4_
//purpose :
//=======================================================================
int mmtrpj4_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *epsi3d,
doublereal *crvlgd,
doublereal *ycvmax,
doublereal *epstrc,
integer *ncfnew)
{
/* Initialized data */
static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
1.05299572648705464724876659688996,
1.0949715351434178709281698645813,
1.15078388379719068145021100764647,
1.2094863084718701596278219811869,
1.26806623151369531323304177532868,
1.32549784426476978866302826176202,
1.38142537365039019558329304432581,
1.43575531950773585146867625840552,
1.48850442653629641402403231015299,
1.53973611681876234549146350844736,
1.58953193485272191557448229046492,
1.63797820416306624705258190017418,
1.68515974143594899185621942934906,
1.73115699602477936547107755854868,
1.77604489805513552087086912113251,
1.81989256661534438347398400420601,
1.86276344480103110090865609776681,
1.90471563564740808542244678597105,
1.94580231994751044968731427898046,
1.98607219357764450634552790950067,
2.02556989246317857340333585562678,
2.06433638992049685189059517340452,
2.10240936014742726236706004607473,
2.13982350649113222745523925190532,
2.17661085564771614285379929798896,
2.21280102016879766322589373557048,
2.2484214321456956597803794333791,
2.28349755104077956674135810027654,
2.31805304852593774867640120860446,
2.35210997297725685169643559615022,
2.38568889602346315560143377261814,
2.41880904328694215730192284109322,
2.45148841120796359750021227795539,
2.48374387161372199992570528025315,
2.5155912654873773953959098501893,
2.54704548720896557684101746505398,
2.57812056037881628390134077704127,
2.60882970619319538196517982945269,
2.63918540521920497868347679257107,
2.66919945330942891495458446613851,
2.69888301230439621709803756505788,
2.72824665609081486737132853370048,
2.75730041251405791603760003778285,
2.78605380158311346185098508516203,
2.81451587035387403267676338931454,
2.84269522483114290814009184272637,
2.87060005919012917988363332454033,
2.89823818258367657739520912946934,
2.92561704377132528239806135133273,
2.95274375377994262301217318010209,
2.97962510678256471794289060402033,
3.00626759936182712291041810228171,
3.03267744830655121818899164295959,
3.05886060707437081434964933864149 };
/* System generated locals */
integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer ncut, i__;
doublereal bidon, error;
integer ia, nd;
doublereal bid, eps1;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Lowers the degree of a curve defined on (-1,1) in the direction of */
/* Legendre with a given precision. */
/* KEYWORDS : */
/* ----------- */
/* LEGENDRE, POLYGON, TRONCATION, CURVE, SMOOTHING. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
/* NDIMEN : Dimension of the space. */
/* NCOEFF : Degree +1 of the polynom. */
/* EPSI3D : Precision required for the approximation. */
/* CRVLGD : The curve which wishes to lower the degree. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* YCVMAX : Auxiliary table (max error on each dimension).
*/
/* EPSTRC : Precision of the approximation. */
/* NCFNEW : Degree +1 of the resulting polynom. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* Parameter adjustments */
--ycvmax;
crvlgd_dim1 = *ncofmx;
crvlgd_offset = crvlgd_dim1 + 1;
crvlgd -= crvlgd_offset;
/* Function Body */
/* Minimum degree that can be reached : Stop at IA (RBD). -------------
*/
ia = 4;
*ncfnew = ia;
/* Init for error calculation. */
i__1 = *ndimen;
for (i__ = 1; i__ <= i__1; ++i__) {
ycvmax[i__] = 0.;
/* L100: */
}
*epstrc = 0.;
error = 0.;
/* Cutting of coefficients. */
ncut = ia + 1;
/* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
*/
i__1 = ncut;
for (i__ = *ncoeff; i__ >= i__1; --i__) {
/* Factor of renormalization. */
bidon = xmaxj[i__ - ncut];
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
bidon;
/* L310: */
}
/* Stop cutting if the norm becomes too great. */
error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
if (error > *epsi3d) {
*ncfnew = i__;
goto L400;
}
/* -- Error max cumulated when the I-eme coeff is removed. */
*epstrc = error;
/* L300: */
}
/* ------- Cutting of zero coeffs of the pole of interpolation (RBD) -------
*/
L400:
if (*ncfnew == ia) {
AdvApp2Var_MathBase::mmeps1_(&eps1);
for (i__ = ia; i__ >= 2; --i__) {
bid = 0.;
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
/* L600: */
}
if (bid > eps1) {
*ncfnew = i__;
goto L9999;
}
/* L500: */
}
/* --- If all coeffs can be removed, this is a point. */
*ncfnew = 1;
}
/* --------------------------------- End --------------------------------
*/
L9999:
return 0;
} /* mmtrpj4_ */
//=======================================================================
//function : mmtrpj6_
//purpose :
//=======================================================================
int mmtrpj6_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *epsi3d,
doublereal *crvlgd,
doublereal *ycvmax,
doublereal *epstrc,
integer *ncfnew)
{
/* Initialized data */
static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
1.11626917091567929907256116528817,
1.1327140810290884106278510474203,
1.1679452722668028753522098022171,
1.20910611986279066645602153641334,
1.25228283758701572089625983127043,
1.29591971597287895911380446311508,
1.3393138157481884258308028584917,
1.3821288728999671920677617491385,
1.42420414683357356104823573391816,
1.46546895108549501306970087318319,
1.50590085198398789708599726315869,
1.54550385142820987194251585145013,
1.58429644271680300005206185490937,
1.62230484071440103826322971668038,
1.65955905239130512405565733793667,
1.69609056468292429853775667485212,
1.73193098017228915881592458573809,
1.7671112206990325429863426635397,
1.80166107681586964987277458875667,
1.83560897003644959204940535551721,
1.86898184653271388435058371983316,
1.90180515174518670797686768515502,
1.93410285411785808749237200054739,
1.96589749778987993293150856865539,
1.99721027139062501070081653790635,
2.02806108474738744005306947877164,
2.05846864831762572089033752595401,
2.08845055210580131460156962214748,
2.11802334209486194329576724042253,
2.14720259305166593214642386780469,
2.17600297710595096918495785742803,
2.20443832785205516555772788192013,
2.2325216999457379530416998244706,
2.2602654243075083168599953074345,
2.28768115912702794202525264301585,
2.3147799369092684021274946755348,
2.34157220782483457076721300512406,
2.36806787963276257263034969490066,
2.39427635443992520016789041085844,
2.42020656255081863955040620243062,
2.44586699364757383088888037359254,
2.47126572552427660024678584642791,
2.49641045058324178349347438430311,
2.52130850028451113942299097584818,
2.54596686772399937214920135190177,
2.5703922285006754089328998222275,
2.59459096001908861492582631591134,
2.61856915936049852435394597597773,
2.64233265984385295286445444361827,
2.66588704638685848486056711408168,
2.68923766976735295746679957665724,
2.71238965987606292679677228666411 };
/* System generated locals */
integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
doublereal d__1;
/* Local variables */
integer ncut, i__;
doublereal bidon, error;
integer ia, nd;
doublereal bid, eps1;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Lowers the degree of a curve defined on (-1,1) in the direction of */
/* Legendre to a given precision. */
/* KEYWORDS : */
/* ----------- */
/* LEGENDRE,POLYGON,TRUNCATION,CURVE,SMOOTHING. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
/* NDIMEN : Dimension of the space. */
/* NCOEFF : Degree +1 of the polynom. */
/* EPSI3D : Precision required for the approximation. */
/* CRVLGD : The curve the degree which of will be lowered. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* YCVMAX : Auxiliary table (max error on each dimension). */
/* EPSTRC : Precision of the approximation. */
/* NCFNEW : Degree +1 of the resulting polynom. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* Parameter adjustments */
--ycvmax;
crvlgd_dim1 = *ncofmx;
crvlgd_offset = crvlgd_dim1 + 1;
crvlgd -= crvlgd_offset;
/* Function Body */
/* Minimum degree that can be reached : Stop at IA (RBD). -------------
*/
ia = 6;
*ncfnew = ia;
/* Init for error calculation. */
i__1 = *ndimen;
for (i__ = 1; i__ <= i__1; ++i__) {
ycvmax[i__] = 0.;
/* L100: */
}
*epstrc = 0.;
error = 0.;
/* Cutting of coefficients. */
ncut = ia + 1;
/* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
*/
i__1 = ncut;
for (i__ = *ncoeff; i__ >= i__1; --i__) {
/* Factor of renormalization. */
bidon = xmaxj[i__ - ncut];
i__2 = *ndimen;
for (nd = 1; nd <= i__2; ++nd) {
ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
bidon;
/* L310: */
}
/* Stop cutting if the norm becomes too great. */
error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
if (error > *epsi3d) {
*ncfnew = i__;
goto L400;
}
/* --- Max error cumulated when the I-th coeff is removed. */
*epstrc = error;
/* L300: */
}
/* ------- Cutting of zero coeff. of the pole of interpolation (RBD) -------
*/
L400:
if (*ncfnew == ia) {
AdvApp2Var_MathBase::mmeps1_(&eps1);
for (i__ = ia; i__ >= 2; --i__) {
bid = 0.;
i__1 = *ndimen;
for (nd = 1; nd <= i__1; ++nd) {
bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
/* L600: */
}
if (bid > eps1) {
*ncfnew = i__;
goto L9999;
}
/* L500: */
}
/* --- If all coeffs can be removed, this is a point. */
*ncfnew = 1;
}
/* --------------------------------- End --------------------------------
*/
L9999:
return 0;
} /* mmtrpj6_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmtrpjj_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmtrpjj_(integer *ncofmx,
integer *ndimen,
integer *ncoeff,
doublereal *epsi3d,
integer *iordre,
doublereal *crvlgd,
doublereal *ycvmax,
doublereal *errmax,
integer *ncfnew)
{
/* System generated locals */
integer crvlgd_dim1, crvlgd_offset;
/* Local variables */
integer ia;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Lower the degree of a curve defined on (-1,1) in the direction of */
/* Legendre with a given precision. */
/* KEYWORDS : */
/* ----------- */
/* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOFMX : Max Nb coeff. of the curve (dimensioning). */
/* NDIMEN : Dimension of the space. */
/* NCOEFF : Degree +1 of the polynom. */
/* EPSI3D : Precision required for the approximation. */
/* IORDRE : Order of continuity at the extremities. */
/* CRVLGD : The curve the degree which of should be lowered. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* ERRMAX : Precision of the approximation. */
/* NCFNEW : Degree +1 of the resulting polynom. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* Parameter adjustments */
--ycvmax;
crvlgd_dim1 = *ncofmx;
crvlgd_offset = crvlgd_dim1 + 1;
crvlgd -= crvlgd_offset;
/* Function Body */
ia = (*iordre + 1) << 1;
if (ia == 0) {
mmtrpj0_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
ycvmax[1], errmax, ncfnew);
} else if (ia == 2) {
mmtrpj2_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
ycvmax[1], errmax, ncfnew);
} else if (ia == 4) {
mmtrpj4_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
ycvmax[1], errmax, ncfnew);
} else {
mmtrpj6_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
ycvmax[1], errmax, ncfnew);
}
/* ------------------------ End -----------------------------------------
*/
return 0;
} /* mmtrpjj_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmunivt_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmunivt_(integer *ndimen,
doublereal *vector,
doublereal *vecnrm,
doublereal *epsiln,
integer *iercod)
{
doublereal c_b2 = 10.;
/* System generated locals */
integer i__1;
doublereal d__1;
/* Local variables */
integer nchif, iunit = 1, izero;
doublereal vnorm;
integer ii;
doublereal bid;
doublereal eps0;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* CALCULATE THE NORMAL VECTOR BASING ON ANY VECTOR */
/* WITH PRECISION GIVEN BY THE USER. */
/* KEYWORDS : */
/* ----------- */
/* ALL, MATH_ACCES :: */
/* VECTEUR&, NORMALISATION, &VECTEUR */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMEN : DIMENSION OF THE SPACE */
/* VECTOR : VECTOR TO BE NORMED */
/* EPSILN : EPSILON BELOW WHICH IT IS CONSIDERED THAT THE */
/* NORM OF THE VECTOR IS NULL. IF EPSILN<=0, A DEFAULT VALUE */
/* IS IMPOSED (10.D-17 ON VAX). */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* VECNRM : NORMED VECTOR */
/* IERCOD 101 : THE VECTOR IS NULL UP TO EPSILN. */
/* 0 : OK. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* VECTOR and VECNRM can be identic. */
/* The norm of vector is calculated and each component is divided by */
/* this norm. After this it is checked if all componentes of the */
/* vector except for one cost 0 with machine precision. In */
/* this case the quasi-null components are set to 0.D0. */
/* > */
/* ***********************************************************************
*/
/* Parameter adjustments */
--vecnrm;
--vector;
/* Function Body */
*iercod = 0;
/* -------- Precision by default : zero machine 10.D-17 on Vax ------
*/
AdvApp2Var_SysBase::maovsr8_(&nchif);
if (*epsiln <= 0.) {
i__1 = -nchif;
eps0 = AdvApp2Var_MathBase::pow__di(&c_b2, &i__1);
} else {
eps0 = *epsiln;
}
/* ------------------------- Calculation of the norm --------------------
*/
vnorm = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vector[1]);
if (vnorm <= eps0) {
AdvApp2Var_SysBase::mvriraz_(ndimen, &vecnrm[1]);
*iercod = 101;
goto L9999;
}
/* ---------------------- Calculation of the vector norm ---------------
*/
izero = 0;
i__1 = (-nchif - 1) / 2;
eps0 = AdvApp2Var_MathBase::pow__di(&c_b2, &i__1);
i__1 = *ndimen;
for (ii = 1; ii <= i__1; ++ii) {
vecnrm[ii] = vector[ii] / vnorm;
if ((d__1 = vecnrm[ii], advapp_abs(d__1)) <= eps0) {
++izero;
} else {
iunit = ii;
}
/* L20: */
}
/* ------ Case when all coordinates except for one are almost null ----
*/
/* ------------- then one of coordinates costs 1.D0 or -1.D0 --------
*/
if (izero == *ndimen - 1) {
bid = vecnrm[iunit];
i__1 = *ndimen;
for (ii = 1; ii <= i__1; ++ii) {
vecnrm[ii] = 0.;
/* L30: */
}
if (bid > 0.) {
vecnrm[iunit] = 1.;
} else {
vecnrm[iunit] = -1.;
}
}
/* -------------------------------- The end -----------------------------
*/
L9999:
return 0;
} /* mmunivt_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmveps3_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmveps3_(doublereal *eps03)
{
/* Initialized data */
static char nomprg[8+1] = "MMEPS1 ";
integer ibb;
/************************************************************************
*******/
/* FUNCTION : */
/* ---------- */
/* Extraction of EPS1 from COMMON MPRCSN. */
/* KEYWORDS : */
/* ----------- */
/* MPRCSN,PRECISON,EPS3. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* Humm. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* EPS3 : space zero of the denominator (10**-9) */
/* EPS3 should value 10**-15 */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* GIVES TOLERANCES OF NULLITY IN STRIM */
/* AND LIMITS OF ITERATIVE PROCESSES */
/* GENERAL CONTEXT, MODIFIABLE BY THE UTILISER */
/* KEYWORDS : */
/* ----------- */
/* PARAMETER , TOLERANCE */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* INITIALISATION : PROFILE , **VIA MPRFTX** AT INPUT IN STRIM*/
/* LOADING OF DEFAULT VALUES OF THE PROFILE IN MPRFTX AT INPUT*/
/* IN STRIM. THEY ARE PRESERVED IN THE LOCAL VARIABLES OF MPRFTX */
/* RESET DEFAULT VALUES : MDFINT */
/* MODIFICATION INTERACTIVE BY THE USER : MDBINT */
/* ACCESS FUNCTION : MMEPS1 ... EPS1 */
/* MEPSPB ... EPS3,EPS4 */
/* MEPSLN ... EPS2, NITERM , NITERR */
/* MEPSNR ... EPS2 , NITERM */
/* MITERR ... NITERR */
/* > */
/* ***********************************************************************
*/
/* NITERM : MAX NB OF ITERATIONS */
/* NITERR : NB OF RAPID ITERATIONS */
/* EPS1 : TOLERANCE OF 3D NULL DISTANCE */
/* EPS2 : TOLERANCE OF ZERO PARAMETRIC DISTANCE */
/* EPS3 : TOLERANCE TO AVOID DIVISION BY 0.. */
/* EPS4 : TOLERANCE ANGULAR */
/* ***********************************************************************
*/
ibb = AdvApp2Var_SysBase::mnfndeb_();
if (ibb >= 5) {
AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
}
*eps03 = mmprcsn_.eps3;
return 0;
} /* mmveps3_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmvncol_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mmvncol_(integer *ndimen,
doublereal *vecin,
doublereal *vecout,
integer *iercod)
{
/* System generated locals */
integer i__1;
/* Local variables */
logical ldbg;
integer d__;
doublereal vaux1[3], vaux2[3];
logical colin;
doublereal valaux;
integer aux;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* CALCULATE A VECTOR NON-COLINEAR TO A GIVEN NON-NULL VECTOR */
/* KEYWORDS : */
/* ----------- */
/* PUBLIC, VECTOR, FREE */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* ndimen : dimension of the space */
/* vecin : input vector */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* vecout : vector non colinear to vecin */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* INITIALISATIONS */
/* ***********************************************************************
*/
/* Parameter adjustments */
--vecout;
--vecin;
/* Function Body */
ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
if (ldbg) {
AdvApp2Var_SysBase::mgenmsg_("MMVNCOL", 7L);
}
*iercod = 0;
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
if (*ndimen <= 1 || *ndimen > 3) {
goto L9101;
}
d__ = 1;
aux = 0;
while(d__ <= *ndimen) {
if (vecin[d__] == 0.) {
++aux;
}
++d__;
}
if (aux == *ndimen) {
goto L9101;
}
for (d__ = 1; d__ <= 3; ++d__) {
vaux1[d__ - 1] = 0.;
}
i__1 = *ndimen;
for (d__ = 1; d__ <= i__1; ++d__) {
vaux1[d__ - 1] = vecin[d__];
vaux2[d__ - 1] = vecin[d__];
}
colin = TRUE_;
d__ = 0;
while(colin) {
++d__;
if (d__ > 3) {
goto L9101;
}
vaux2[d__ - 1] += 1;
valaux = vaux1[1] * vaux2[2] - vaux1[2] * vaux2[1];
if (valaux == 0.) {
valaux = vaux1[2] * vaux2[0] - vaux1[0] * vaux2[2];
if (valaux == 0.) {
valaux = vaux1[0] * vaux2[1] - vaux1[1] * vaux2[0];
if (valaux != 0.) {
colin = FALSE_;
}
} else {
colin = FALSE_;
}
} else {
colin = FALSE_;
}
}
if (colin) {
goto L9101;
}
i__1 = *ndimen;
for (d__ = 1; d__ <= i__1; ++d__) {
vecout[d__] = vaux2[d__ - 1];
}
goto L9999;
/* ***********************************************************************
*/
/* ERROR PROCESSING */
/* ***********************************************************************
*/
L9101:
*iercod = 1;
goto L9999;
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
L9999:
AdvApp2Var_SysBase::maermsg_("MMVNCOL", iercod, 7L);
if (ldbg) {
AdvApp2Var_SysBase::mgsomsg_("MMVNCOL", 7L);
}
return 0 ;
} /* mmvncol_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mmwprcs_
//purpose :
//=======================================================================
void AdvApp2Var_MathBase::mmwprcs_(doublereal *epsil1,
doublereal *epsil2,
doublereal *epsil3,
doublereal *epsil4,
integer *niter1,
integer *niter2)
{
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* ACCESS IN WRITING FOR COMMON MPRCSN */
/* KEYWORDS : */
/* ----------- */
/* WRITING */
/* INPUT ARGUMENTS : */
/* -------------------- */
/* EPSIL1 : TOLERANCE OF 3D NULL DISTANCE */
/* EPSIL2 : TOLERANCE OF PARAMETRIC NULL DISTANCE */
/* EPSIL3 : TOLERANCE TO AVOID DIVISION BY 0.. */
/* EPSIL4 : ANGULAR TOLERANCE */
/* NITER1 : MAX NB OF ITERATIONS */
/* NITER2 : NB OF RAPID ITERATIONS */
/* OUTPUT ARGUMENTS : */
/* --------------------- */
/* NONE */
/* COMMONS USED : */
/* ------------------ */
/* REFERENCES CALLED : */
/* --------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* INITIALIZATIONS */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* GIVES TOLERANCES OF NULLITY IN STRIM */
/* AND LIMITS OF ITERATIVE PROCESSES */
/* GENERAL CONTEXT, MODIFIABLE BY THE UTILISER */
/* KEYWORDS : */
/* ----------- */
/* PARAMETER , TOLERANCE */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* INITIALISATION : PROFILE , **VIA MPRFTX** AT INPUT IN STRIM*/
/* LOADING OF DEFAULT VALUES OF THE PROFILE IN MPRFTX AT INPUT*/
/* IN STRIM. THEY ARE PRESERVED IN THE LOCAL VARIABLES OF MPRFTX */
/* RESET DEFAULT VALUES : MDFINT */
/* MODIFICATION INTERACTIVE BY THE USER : MDBINT */
/* ACCESS FUNCTION : MMEPS1 ... EPS1 */
/* MEPSPB ... EPS3,EPS4 */
/* MEPSLN ... EPS2, NITERM , NITERR */
/* MEPSNR ... EPS2 , NITERM */
/* MITERR ... NITERR */
/* > */
/* ***********************************************************************
*/
/* NITERM : MAX NB OF ITERATIONS */
/* NITERR : NB OF RAPID ITERATIONS */
/* EPS1 : TOLERANCE OF 3D NULL DISTANCE */
/* EPS2 : TOLERANCE OF ZERO PARAMETRIC DISTANCE */
/* EPS3 : TOLERANCE TO AVOID DIVISION BY 0.. */
/* EPS4 : TOLERANCE ANGULAR */
/* ***********************************************************************
*/
mmprcsn_.eps1 = *epsil1;
mmprcsn_.eps2 = *epsil2;
mmprcsn_.eps3 = *epsil3;
mmprcsn_.eps4 = *epsil4;
mmprcsn_.niterm = *niter1;
mmprcsn_.niterr = *niter2;
return ;
} /* mmwprcs_ */
//=======================================================================
//function : AdvApp2Var_MathBase::pow__di
//purpose :
//=======================================================================
doublereal AdvApp2Var_MathBase::pow__di (doublereal *x,
integer *n)
{
register integer ii ;
doublereal result ;
integer absolute ;
result = 1.0e0 ;
if ( *n > 0 ) {absolute = *n;}
else {absolute = -*n;}
/* System generated locals */
for(ii = 0 ; ii < absolute ; ii++) {
result *= *x ;
}
if (*n < 0) {
result = 1.0e0 / result ;
}
return result ;
}
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Calculate integer function power not obligatory in the most efficient way ;
*/
/* KEYWORDS : */
/* ----------- */
/* POWER */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* X : argument of X**N */
/* N : power */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* return X**N */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************/
//=======================================================================
//function : pow__ii
//purpose :
//=======================================================================
integer pow__ii(integer *x,
integer *n)
{
register integer ii ;
integer result ;
integer absolute ;
result = 1 ;
if ( *n > 0 ) {absolute = *n;}
else {absolute = -*n;}
/* System generated locals */
for(ii = 0 ; ii < absolute ; ii++) {
result *= *x ;
}
if (*n < 0) {
result = 1 / result ;
}
return result ;
}
/* **********************************************************************
*/
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Calculate integer function power not obligatory in the most efficient way ;
*/
/* KEYWORDS : */
/* ----------- */
/* POWER */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* X : argument of X**N */
/* N : power */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* return X**N */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************/
//=======================================================================
//function : AdvApp2Var_MathBase::msc_
//purpose :
//=======================================================================
doublereal AdvApp2Var_MathBase::msc_(integer *ndimen,
doublereal *vecte1,
doublereal *vecte2)
{
/* System generated locals */
integer i__1;
doublereal ret_val;
/* Local variables */
integer i__;
doublereal x;
/************************************************************************
*******/
/* FUNCTION : */
/* ---------- */
/* Calculate the scalar product of 2 vectors in the space */
/* of dimension NDIMEN. */
/* KEYWORDS : */
/* ----------- */
/* PRODUCT MSCALAIRE. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMEN : Dimension of the space. */
/* VECTE1,VECTE2: Vectors. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ----------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* > */
/* ***********************************************************************
*/
/* PRODUIT MSCALAIRE */
/* Parameter adjustments */
--vecte2;
--vecte1;
/* Function Body */
x = 0.;
i__1 = *ndimen;
for (i__ = 1; i__ <= i__1; ++i__) {
x += vecte1[i__] * vecte2[i__];
/* L100: */
}
ret_val = x;
/* ----------------------------------- THE END --------------------------
*/
return ret_val;
} /* msc_ */
//=======================================================================
//function : mvcvin2_
//purpose :
//=======================================================================
int mvcvin2_(integer *ncoeff,
doublereal *crvold,
doublereal *crvnew,
integer *iercod)
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
integer m1jm1, ncfm1, j, k;
doublereal bid;
doublereal cij1, cij2;
/************************************************************************
*******/
/* FONCTION : */
/* ---------- */
/* INVERSION OF THE PARAMETERS ON CURVE 2D. */
/* KEYWORDS : */
/* ----------- */
/* CURVE,2D,INVERSION,PARAMETER. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOEFF : NB OF COEFF OF THE CURVE. */
/* CRVOLD : CURVE OF ORIGIN */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* CRVNEW : THE RESULTING CURVE AFTER CHANGE OF T BY 1-T */
/* IERCOD : 0 OK, */
/* 10 NB OF COEFF NULL OR TOO GREAT. */
/* COMMONS USED : */
/* ---------------- */
/* MCCNP */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Neant */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* THE FOLLOWING CALL IS ABSOLUTELY LEGAL : */
/* CALL MVCVIN2(NCOEFF,CURVE,CURVE,IERCOD), THE TABLE CURVE */
/* BECOMES INPUT AND OUTPUT ARGUMENT (RBD). */
/* BECAUSE OF MCCNP, THE NB OF COEFF OF THE CURVE IS LIMITED TO */
/* NDGCNP+1 = 61. */
/* > */
/* ***********************************************************************
*/
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Serves to provide coefficients of the binome (triangle of Pascal). */
/* KEYWORDS : */
/* ----------- */
/* Coeff of binome from 0 to 60. read only . init par block data */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* The coefficients of the binome form a triangular matrix. */
/* This matrix is completed in table CNP by transposition. */
/* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
/* Initialization is done by block-data MMLLL09.RES, */
/* created by program MQINICNP.FOR (see the team (AC) ). */
/* > */
/* **********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
crvnew -= 3;
crvold -= 3;
/* Function Body */
if (*ncoeff < 1 || *ncoeff - 1 > 60) {
*iercod = 10;
goto L9999;
}
*iercod = 0;
/* CONSTANT TERM OF THE NEW CURVE */
cij1 = crvold[3];
cij2 = crvold[4];
i__1 = *ncoeff;
for (k = 2; k <= i__1; ++k) {
cij1 += crvold[(k << 1) + 1];
cij2 += crvold[(k << 1) + 2];
}
crvnew[3] = cij1;
crvnew[4] = cij2;
if (*ncoeff == 1) {
goto L9999;
}
/* INTERMEDIARY POWERS OF THE PARAMETER */
ncfm1 = *ncoeff - 1;
m1jm1 = 1;
i__1 = ncfm1;
for (j = 2; j <= i__1; ++j) {
m1jm1 = -m1jm1;
cij1 = crvold[(j << 1) + 1];
cij2 = crvold[(j << 1) + 2];
i__2 = *ncoeff;
for (k = j + 1; k <= i__2; ++k) {
bid = mmcmcnp_.cnp[k - 1 + (j - 1) * 61];
cij1 += crvold[(k << 1) + 1] * bid;
cij2 += crvold[(k << 1) + 2] * bid;
}
crvnew[(j << 1) + 1] = cij1 * m1jm1;
crvnew[(j << 1) + 2] = cij2 * m1jm1;
}
/* TERM OF THE HIGHEST DEGREE */
crvnew[(*ncoeff << 1) + 1] = -crvold[(*ncoeff << 1) + 1] * m1jm1;
crvnew[(*ncoeff << 1) + 2] = -crvold[(*ncoeff << 1) + 2] * m1jm1;
L9999:
if (*iercod > 0) {
AdvApp2Var_SysBase::maermsg_("MVCVIN2", iercod, 7L);
}
return 0 ;
} /* mvcvin2_ */
//=======================================================================
//function : mvcvinv_
//purpose :
//=======================================================================
int mvcvinv_(integer *ncoeff,
doublereal *crvold,
doublereal *crvnew,
integer *iercod)
{
/* System generated locals */
integer i__1, i__2;
/* Local variables */
integer m1jm1, ncfm1, j, k;
doublereal bid;
//extern /* Subroutine */ int maermsg_();
doublereal cij1, cij2, cij3;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* INVERSION OF THE PARAMETER ON A CURBE 3D (I.E. INVERSION */
/* OF THE DIRECTION OF PARSING). */
/* KEYWORDS : */
/* ----------- */
/* CURVE,INVERSION,PARAMETER. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOEFF : NB OF COEFF OF THE CURVE. */
/* CRVOLD : CURVE OF ORIGIN */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* CRVNEW : RESULTING CURVE AFTER CHANGE OF T INTO 1-T */
/* IERCOD : 0 OK, */
/* 10 NB OF COEFF NULL OR TOO GREAT. */
/* COMMONS USED : */
/* ---------------- */
/* MCCNP */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Neant */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* THE FOLLOWING CALL IS ABSOLUTELY LEGAL : */
/* CALL MVCVINV(NCOEFF,CURVE,CURVE,IERCOD), TABLE CURVE */
/* BECOMES INPUT AND OUTPUT ARGUMENT (RBD). */
/* THE NUMBER OF COEFF OF THE CURVE IS LIMITED TO NDGCNP+1 = 61 */
/* BECAUSE OF USE OF COMMON MCCNP. */
/* > */
/* ***********************************************************************
*/
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* Serves to provide the binomial coefficients (triangle of Pascal). */
/* KEYWORDS : */
/* ----------- */
/* Binomial Coeff from 0 to 60. read only . init par block data */
/* DEMSCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* The binomial coefficients form a triangular matrix. */
/* This matrix is completed in table CNP by its transposition. */
/* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
/* Initialisation is done by block-data MMLLL09.RES, */
/* created by program MQINICNP.FOR (see the team (AC) ). */
/* > */
/* **********************************************************************
*/
/* ***********************************************************************
*/
/* Parameter adjustments */
crvnew -= 4;
crvold -= 4;
/* Function Body */
if (*ncoeff < 1 || *ncoeff - 1 > 60) {
*iercod = 10;
goto L9999;
}
*iercod = 0;
/* CONSTANT TERM OF THE NEW CURVE */
cij1 = crvold[4];
cij2 = crvold[5];
cij3 = crvold[6];
i__1 = *ncoeff;
for (k = 2; k <= i__1; ++k) {
cij1 += crvold[k * 3 + 1];
cij2 += crvold[k * 3 + 2];
cij3 += crvold[k * 3 + 3];
/* L30: */
}
crvnew[4] = cij1;
crvnew[5] = cij2;
crvnew[6] = cij3;
if (*ncoeff == 1) {
goto L9999;
}
/* INTERMEDIARY POWER OF THE PARAMETER */
ncfm1 = *ncoeff - 1;
m1jm1 = 1;
i__1 = ncfm1;
for (j = 2; j <= i__1; ++j) {
m1jm1 = -m1jm1;
cij1 = crvold[j * 3 + 1];
cij2 = crvold[j * 3 + 2];
cij3 = crvold[j * 3 + 3];
i__2 = *ncoeff;
for (k = j + 1; k <= i__2; ++k) {
bid = mmcmcnp_.cnp[k - 1 + (j - 1) * 61];
cij1 += crvold[k * 3 + 1] * bid;
cij2 += crvold[k * 3 + 2] * bid;
cij3 += crvold[k * 3 + 3] * bid;
/* L40: */
}
crvnew[j * 3 + 1] = cij1 * m1jm1;
crvnew[j * 3 + 2] = cij2 * m1jm1;
crvnew[j * 3 + 3] = cij3 * m1jm1;
/* L50: */
}
/* TERM OF THE HIGHEST DEGREE */
crvnew[*ncoeff * 3 + 1] = -crvold[*ncoeff * 3 + 1] * m1jm1;
crvnew[*ncoeff * 3 + 2] = -crvold[*ncoeff * 3 + 2] * m1jm1;
crvnew[*ncoeff * 3 + 3] = -crvold[*ncoeff * 3 + 3] * m1jm1;
L9999:
AdvApp2Var_SysBase::maermsg_("MVCVINV", iercod, 7L);
return 0;
} /* mvcvinv_ */
//=======================================================================
//function : mvgaus0_
//purpose :
//=======================================================================
int mvgaus0_(integer *kindic,
doublereal *urootl,
doublereal *hiltab,
integer *nbrval,
integer *iercod)
{
/* System generated locals */
integer i__1;
/* Local variables */
doublereal tamp[40];
integer ndegl, kg, ii;
/* **********************************************************************
*/
/* FUNCTION : */
/* -------- */
/* Loading of a degree gives roots of LEGENDRE polynom */
/* DEFINED on [-1,1] and weights of Gauss quadrature formulas */
/* (based on corresponding LAGRANGIAN interpolators). */
/* The symmetry relative to 0 is used between [-1,0] and [0,1]. */
/* KEYWORDS : */
/* --------- */
/* . VOLUMIC, LEGENDRE, LAGRANGE, GAUSS */
/* INPUT ARGUMENTSE : */
/* ------------------ */
/* KINDIC : Takes values from 1 to 10 depending of the degree */
/* of the used polynom. */
/* The degree of the polynom is equal to 4 k, i.e. 4, 8, */
/* 12, 16, 20, 24, 28, 32, 36 and 40. */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* UROOTL : Roots of LEGENDRE polynom in domain [1,0] */
/* given in decreasing order. For domain [-1,0], it is */
/* necessary to take the opposite values. */
/* HILTAB : LAGRANGE interpolators associated to roots. For */
/* opposed roots, interpolatorsare equal. */
/* NBRVAL : Nb of coefficients. Is equal to the half of degree */
/* depending on the symmetry (i.e. 2*KINDIC). */
/* IERCOD : Error code: */
/* < 0 ==> Attention - Warning */
/* =-1 ==> Value of false KINDIC. NBRVAL is forced to 20 */
/* (order 40) */
/* = 0 ==> Everything is OK */
/* COMMON USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* --------------------------------- */
/* If KINDIC is not correct (i.e < 1 or > 10), the degree is set */
/* to 40 directly (ATTENTION to overload - to avoid it, */
/* preview UROOTL and HILTAB dimensioned at least to 20). */
/* The value of coefficients was calculated with quadruple precision */
/* by JJM with help of GD. */
/* Checking of roots was done by GD. */
/* See detailed explications on the listing */
/* > */
/* **********************************************************************
*/
/* ------------------------------------ */
/* ****** Test validity of KINDIC ** */
/* ------------------------------------ */
/* Parameter adjustments */
--hiltab;
--urootl;
/* Function Body */
*iercod = 0;
kg = *kindic;
if (kg < 1 || kg > 10) {
kg = 10;
*iercod = -1;
}
*nbrval = kg << 1;
ndegl = *nbrval << 1;
/* ----------------------------------------------------------------------
*/
/* ****** Load NBRVAL positive roots depending on the degree **
*/
/* ----------------------------------------------------------------------
*/
/* ATTENTION : Sign minus (-) in the loop is intentional. */
mmextrl_(&ndegl, tamp);
i__1 = *nbrval;
for (ii = 1; ii <= i__1; ++ii) {
urootl[ii] = -tamp[ii - 1];
/* L100: */
}
/* ------------------------------------------------------------------- */
/* ****** Loading of NBRVAL Gauss weight depending on the degree ** */
/* ------------------------------------------------------------------- */
mmexthi_(&ndegl, tamp);
i__1 = *nbrval;
for (ii = 1; ii <= i__1; ++ii) {
hiltab[ii] = tamp[ii - 1];
/* L200: */
}
/* ------------------------------- */
/* ****** End of sub-program ** */
/* ------------------------------- */
return 0;
} /* mvgaus0_ */
//=======================================================================
//function : mvpscr2_
//purpose :
//=======================================================================
int mvpscr2_(integer *ncoeff,
doublereal *curve2,
doublereal *tparam,
doublereal *pntcrb)
{
/* System generated locals */
integer i__1;
/* Local variables */
integer ndeg, kk;
doublereal xxx, yyy;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* POSITIONING ON CURVE (NCF,2) IN SPACE OF DIMENSION 2. */
/* KEYWORDS : */
/* ----------- */
/* TOUS,MATH_ACCES:: COURBE&,POSITIONNEMENT,&POINT. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOEFF : NUMBER OF COEFFICIENTS OF THE CURVE */
/* CURVE2 : EQUATION OF CURVE 2D */
/* TPARAM : VALUE OF PARAMETER AT GIVEN POINT */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* PNTCRB : COORDINATES OF POINT CORRESPONDING TO PARAMETER */
/* TPARAM ON CURVE 2D CURVE2. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ---------------------- */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* MSCHEMA OF HORNER. */
/* > */
/* **********************************************************************
*/
/* -------- INITIALIZATIONS AND PROCESSING OF PARTICULAR CASES ----------
*/
/* ---> Cas when NCOEFF > 1 (case STANDARD). */
/* Parameter adjustments */
--pntcrb;
curve2 -= 3;
/* Function Body */
if (*ncoeff >= 2) {
goto L1000;
}
/* ---> Case when NCOEFF <= 1. */
if (*ncoeff <= 0) {
pntcrb[1] = 0.;
pntcrb[2] = 0.;
goto L9999;
} else if (*ncoeff == 1) {
pntcrb[1] = curve2[3];
pntcrb[2] = curve2[4];
goto L9999;
}
/* -------------------- MSCHEMA OF HORNER (PARTICULAR CASE) --------------
*/
L1000:
if (*tparam == 1.) {
xxx = 0.;
yyy = 0.;
i__1 = *ncoeff;
for (kk = 1; kk <= i__1; ++kk) {
xxx += curve2[(kk << 1) + 1];
yyy += curve2[(kk << 1) + 2];
/* L100: */
}
goto L5000;
} else if (*tparam == 0.) {
pntcrb[1] = curve2[3];
pntcrb[2] = curve2[4];
goto L9999;
}
/* ---------------------------- MSCHEMA OF HORNER ------------------------
*/
/* ---> TPARAM is different from 1.D0 and 0.D0. */
ndeg = *ncoeff - 1;
xxx = curve2[(*ncoeff << 1) + 1];
yyy = curve2[(*ncoeff << 1) + 2];
for (kk = ndeg; kk >= 1; --kk) {
xxx = xxx * *tparam + curve2[(kk << 1) + 1];
yyy = yyy * *tparam + curve2[(kk << 1) + 2];
/* L200: */
}
goto L5000;
/* ------------------------ RECOVER THE CALCULATED POINT ---------------
*/
L5000:
pntcrb[1] = xxx;
pntcrb[2] = yyy;
/* ------------------------------ THE END -------------------------------
*/
L9999:
return 0;
} /* mvpscr2_ */
//=======================================================================
//function : mvpscr3_
//purpose :
//=======================================================================
int mvpscr3_(integer *ncoeff,
doublereal *curve3,
doublereal *tparam,
doublereal *pntcrb)
{
/* System generated locals */
integer i__1;
/* Local variables */
integer ndeg, kk;
doublereal xxx, yyy, zzz;
/* **********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* POSITIONING ON A CURVE (3,NCF) IN THE SPACE OF DIMENSION 3. */
/* KEYWORDS : */
/* ----------- */
/* TOUS, MATH_ACCES:: COURBE&,POSITIONNEMENT,&POINT. */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NCOEFF : NB OF COEFFICIENTS OF THE CURVE */
/* CURVE3 : EQUATION OF CURVE 3D */
/* TPARAM : VALUE OF THE PARAMETER AT THE GIVEN POINT */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* PNTCRB : COORDINATES OF THE POINT CORRESPONDING TO PARAMETER */
/* TPARAM ON CURVE 3D CURVE3. */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Neant */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* MSCHEMA OF HORNER. */
/* > */
/* **********************************************************************
*/
/* DECLARATIONS */
/* **********************************************************************
*/
/* -------- INITIALISATIONS AND PROCESSING OF PARTICULAR CASES ----------
*/
/* ---> Case when NCOEFF > 1 (cas STANDARD). */
/* Parameter adjustments */
--pntcrb;
curve3 -= 4;
/* Function Body */
if (*ncoeff >= 2) {
goto L1000;
}
/* ---> Case when NCOEFF <= 1. */
if (*ncoeff <= 0) {
pntcrb[1] = 0.;
pntcrb[2] = 0.;
pntcrb[3] = 0.;
goto L9999;
} else if (*ncoeff == 1) {
pntcrb[1] = curve3[4];
pntcrb[2] = curve3[5];
pntcrb[3] = curve3[6];
goto L9999;
}
/* -------------------- MSCHEMA OF HORNER (PARTICULAR CASE) --------------
*/
L1000:
if (*tparam == 1.) {
xxx = 0.;
yyy = 0.;
zzz = 0.;
i__1 = *ncoeff;
for (kk = 1; kk <= i__1; ++kk) {
xxx += curve3[kk * 3 + 1];
yyy += curve3[kk * 3 + 2];
zzz += curve3[kk * 3 + 3];
/* L100: */
}
goto L5000;
} else if (*tparam == 0.) {
pntcrb[1] = curve3[4];
pntcrb[2] = curve3[5];
pntcrb[3] = curve3[6];
goto L9999;
}
/* ---------------------------- MSCHEMA OF HORNER ------------------------
*/
/* ---> Here TPARAM is different from 1.D0 and 0.D0. */
ndeg = *ncoeff - 1;
xxx = curve3[*ncoeff * 3 + 1];
yyy = curve3[*ncoeff * 3 + 2];
zzz = curve3[*ncoeff * 3 + 3];
for (kk = ndeg; kk >= 1; --kk) {
xxx = xxx * *tparam + curve3[kk * 3 + 1];
yyy = yyy * *tparam + curve3[kk * 3 + 2];
zzz = zzz * *tparam + curve3[kk * 3 + 3];
/* L200: */
}
goto L5000;
/* ------------------------ RETURN THE CALCULATED POINT ------------------
*/
L5000:
pntcrb[1] = xxx;
pntcrb[2] = yyy;
pntcrb[3] = zzz;
/* ------------------------------ THE END -------------------------------
*/
L9999:
return 0;
} /* mvpscr3_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mvsheld_
//purpose :
//=======================================================================
int AdvApp2Var_MathBase::mvsheld_(integer *n,
integer *is,
doublereal *dtab,
integer *icle)
{
/* System generated locals */
integer dtab_dim1, dtab_offset, i__1, i__2;
/* Local variables */
integer incr;
doublereal dsave;
integer i3, i4, i5, incrp1;
/************************************************************************
*******/
/* FUNCTION : */
/* ---------- */
/* PARSING OF COLUMNS OF TABLE OF REAL*8 BY SHELL METHOD*/
/* (IN INCREASING ORDER) */
/* KEYWORDS : */
/* ----------- */
/* POINT-ENTRY, PARSING, SHELL */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* N : NUMBER OF COLUMNS OF THE TABLE */
/* IS : NUMBER OF LINE OF THE TABLE */
/* DTAB : TABLE OF REAL*8 TO BE PARSED */
/* ICLE : POSITION OF THE KEY ON THE COLUMN */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* DTAB : PARSED TABLE */
/* COMMONS USED : */
/* ---------------- */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Neant */
/* DESCRIPTION/NOTES/LIMITATIONS : */
/* ----------------------------------- */
/* CLASSIC SHELL METHOD : PARSING BY SERIES */
/* Declaration DTAB(IS, 1) corresponds to DTAB(IS, *) */
/* > */
/* ***********************************************************************
*/
/* Parameter adjustments */
dtab_dim1 = *is;
dtab_offset = dtab_dim1 + 1;
dtab -= dtab_offset;
/* Function Body */
if (*n <= 1) {
goto L9900;
}
/* ------------------------ */
/* INITIALIZATION OF THE SEQUENCE OF INCREMENTS */
/* FIND THE GREATEST INCREMENT SO THAT INCR < N/9 */
incr = 1;
L1001:
if (incr >= *n / 9) {
goto L1002;
}
/* ----------------------------- */
incr = incr * 3 + 1;
goto L1001;
/* LOOP ON INCREMENTS TILL INCR = 1 */
/* PARSING BY SERIES DISTANT FROM INCR */
L1002:
incrp1 = incr + 1;
/* ----------------- */
i__1 = *n;
for (i3 = incrp1; i3 <= i__1; ++i3) {
/* ---------------------- */
/* SET ELEMENT I3 AT ITS PLACE IN THE SERIES */
i4 = i3 - incr;
L1004:
if (i4 < 1) {
goto L1003;
}
/* ------------------------- */
if (dtab[*icle + i4 * dtab_dim1] <= dtab[*icle + (i4 + incr) *
dtab_dim1]) {
goto L1003;
}
i__2 = *is;
for (i5 = 1; i5 <= i__2; ++i5) {
/* ------------------ */
dsave = dtab[i5 + i4 * dtab_dim1];
dtab[i5 + i4 * dtab_dim1] = dtab[i5 + (i4 + incr) * dtab_dim1];
dtab[i5 + (i4 + incr) * dtab_dim1] = dsave;
}
/* -------- */
i4 -= incr;
goto L1004;
L1003:
;
}
/* -------- */
/* PASSAGE TO THE NEXT INCREMENT */
incr /= 3;
if (incr >= 1) {
goto L1002;
}
L9900:
return 0 ;
} /* mvsheld_ */
//=======================================================================
//function : AdvApp2Var_MathBase::mzsnorm_
//purpose :
//=======================================================================
doublereal AdvApp2Var_MathBase::mzsnorm_(integer *ndimen,
doublereal *vecteu)
{
/* System generated locals */
integer i__1;
doublereal ret_val, d__1, d__2;
/* Local variables */
doublereal xsom;
integer i__, irmax;
/* ***********************************************************************
*/
/* FUNCTION : */
/* ---------- */
/* SERVES to calculate the euclidian norm of a vector : */
/* ____________________________ */
/* Z = V V(1)**2 + V(2)**2 + ... */
/* KEYWORDS : */
/* ----------- */
/* SURMFACIQUE, */
/* INPUT ARGUMENTS : */
/* ------------------ */
/* NDIMEN : Dimension of the vector */
/* VECTEU : vector of dimension NDIMEN */
/* OUTPUT ARGUMENTS : */
/* ------------------- */
/* MZSNORM : Value of the euclidian norm of vector VECTEU */
/* COMMONS USED : */
/* ---------------- */
/* .Neant. */
/* REFERENCES CALLED : */
/* ---------------------- */
/* Type Name */
/* R*8 ABS R*8 SQRT */
/* DESCRIPTION/NOTESS/LIMITATIONS : */
/* ----------------------------------- */
/* To limit the risks of overflow, */
/* the term of the strongest absolute value is factorized : */
/* _______________________ */
/* Z = !V(1)! * V 1 + (V(2)/V(1))**2 + ... */
/* > */
/* ***********************************************************************
*/
/* DECLARATIONS */
/* ***********************************************************************
*/
/* ***********************************************************************
*/
/* PROCESSING */
/* ***********************************************************************
*/
/* ___ Find the strongest absolute value term */
/* Parameter adjustments */
--vecteu;
/* Function Body */
irmax = 1;
i__1 = *ndimen;
for (i__ = 2; i__ <= i__1; ++i__) {
if ((d__1 = vecteu[irmax], advapp_abs(d__1)) < (d__2 = vecteu[i__], advapp_abs(d__2)
)) {
irmax = i__;
}
/* L100: */
}
/* ___ Calculate the norme */
if ((d__1 = vecteu[irmax], advapp_abs(d__1)) < 1.) {
xsom = 0.;
i__1 = *ndimen;
for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing 2nd power */
d__1 = vecteu[i__];
xsom += d__1 * d__1;
/* L200: */
}
ret_val = sqrt(xsom);
} else {
xsom = 0.;
i__1 = *ndimen;
for (i__ = 1; i__ <= i__1; ++i__) {
if (i__ == irmax) {
xsom += 1.;
} else {
/* Computing 2nd power */
d__1 = vecteu[i__] / vecteu[irmax];
xsom += d__1 * d__1;
}
/* L300: */
}
ret_val = (d__1 = vecteu[irmax], advapp_abs(d__1)) * sqrt(xsom);
}
/* ***********************************************************************
*/
/* RETURN CALLING PROGRAM */
/* ***********************************************************************
*/
return ret_val;
} /* mzsnorm_ */
View Git Blame Copy Permalink