0023952: Improving thread-safety of intersections, approximations and other modeling...
[occt.git] / src / AdvApp2Var / AdvApp2Var_MathBase.cxx
CommitLineData
b311480e 1// Copyright (c) 1999-2012 OPEN CASCADE SAS
7fd59977 2//
b311480e 3// The content of this file is subject to the Open CASCADE Technology Public
4// License Version 6.5 (the "License"). You may not use the content of this file
5// except in compliance with the License. Please obtain a copy of the License
6// at http://www.opencascade.org and read it completely before using this file.
7//
8// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
9// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
7fd59977 10//
b311480e 11// The Original Code and all software distributed under the License is
12// distributed on an "AS IS" basis, without warranty of any kind, and the
13// Initial Developer hereby disclaims all such warranties, including without
14// limitation, any warranties of merchantability, fitness for a particular
15// purpose or non-infringement. Please see the License for the specific terms
16// and conditions governing the rights and limitations under the License.
17
18// AdvApp2Var_MathBase.cxx
7fd59977 19#include <math.h>
20#include <AdvApp2Var_SysBase.hxx>
21#include <AdvApp2Var_Data_f2c.hxx>
22#include <AdvApp2Var_MathBase.hxx>
23#include <AdvApp2Var_Data.hxx>
24
25// statics
26static
27int mmchole_(integer *mxcoef,
28 integer *dimens,
29 doublereal *amatri,
30 integer *aposit,
31 integer *posuiv,
32 doublereal *chomat,
33 integer *iercod);
34
35
36
37
38static
39int mmrslss_(integer *mxcoef,
40 integer *dimens,
41 doublereal *smatri,
42 integer *sposit,
43 integer *posuiv,
44 doublereal *mscnmbr,
45 doublereal *soluti,
46 integer *iercod);
47
48static
49int mfac_(doublereal *f,
50 integer *n);
51
52static
53int mmaper0_(integer *ncofmx,
54 integer *ndimen,
55 integer *ncoeff,
56 doublereal *crvlgd,
57 integer *ncfnew,
58 doublereal *ycvmax,
59 doublereal *errmax);
60static
61int mmaper2_(integer *ncofmx,
62 integer *ndimen,
63 integer *ncoeff,
64 doublereal *crvjac,
65 integer *ncfnew,
66 doublereal *ycvmax,
67 doublereal *errmax);
68
69static
70int mmaper4_(integer *ncofmx,
71 integer *ndimen,
72 integer *ncoeff,
73 doublereal *crvjac,
74 integer *ncfnew,
75 doublereal *ycvmax,
76 doublereal *errmax);
77
78static
79int mmaper6_(integer *ncofmx,
80 integer *ndimen,
81 integer *ncoeff,
82 doublereal *crvjac,
83 integer *ncfnew,
84 doublereal *ycvmax,
85 doublereal *errmax);
86
87static
88int mmarc41_(integer *ndimax,
89 integer *ndimen,
90 integer *ncoeff,
91 doublereal *crvold,
92 doublereal *upara0,
93 doublereal *upara1,
94 doublereal *crvnew,
95 integer *iercod);
96
97static
98int mmatvec_(integer *nligne,
99 integer *ncolon,
100 integer *gposit,
101 integer *gnstoc,
102 doublereal *gmatri,
103 doublereal *vecin,
104 integer *deblig,
105 doublereal *vecout,
106 integer *iercod);
107
108static
109int mmcvstd_(integer *ncofmx,
110 integer *ndimax,
111 integer *ncoeff,
112 integer *ndimen,
113 doublereal *crvcan,
114 doublereal *courbe);
115
116static
117int mmdrvcb_(integer *ideriv,
118 integer *ndim,
119 integer *ncoeff,
120 doublereal *courbe,
121 doublereal *tparam,
122 doublereal *tabpnt,
123 integer *iercod);
124
125static
126int mmexthi_(integer *ndegre,
127 doublereal *hwgaus);
128
129static
130int mmextrl_(integer *ndegre,
131 doublereal *rootlg);
132
133
134
135static
136int mmherm0_(doublereal *debfin,
137 integer *iercod);
138
139static
140int mmherm1_(doublereal *debfin,
141 integer *ordrmx,
142 integer *iordre,
143 doublereal *hermit,
144 integer *iercod);
145static
146int mmloncv_(integer *ndimax,
147 integer *ndimen,
148 integer *ncoeff,
149 doublereal *courbe,
150 doublereal *tdebut,
151 doublereal *tfinal,
152 doublereal *xlongc,
153 integer *iercod);
154static
155int mmpojac_(doublereal *tparam,
156 integer *iordre,
157 integer *ncoeff,
158 integer *nderiv,
159 doublereal *valjac,
160 integer *iercod);
161
162static
163int mmrslw_(integer *normax,
164 integer *nordre,
165 integer *ndimen,
166 doublereal *epspiv,
167 doublereal *abmatr,
168 doublereal *xmatri,
169 integer *iercod);
170static
171int mmtmave_(integer *nligne,
172 integer *ncolon,
173 integer *gposit,
174 integer *gnstoc,
175 doublereal *gmatri,
176 doublereal *vecin,
177 doublereal *vecout,
178 integer *iercod);
179static
180int mmtrpj0_(integer *ncofmx,
181 integer *ndimen,
182 integer *ncoeff,
183 doublereal *epsi3d,
184 doublereal *crvlgd,
185 doublereal *ycvmax,
186 doublereal *epstrc,
187 integer *ncfnew);
188static
189int mmtrpj2_(integer *ncofmx,
190 integer *ndimen,
191 integer *ncoeff,
192 doublereal *epsi3d,
193 doublereal *crvlgd,
194 doublereal *ycvmax,
195 doublereal *epstrc,
196 integer *ncfnew);
197
198static
199int mmtrpj4_(integer *ncofmx,
200 integer *ndimen,
201 integer *ncoeff,
202 doublereal *epsi3d,
203 doublereal *crvlgd,
204 doublereal *ycvmax,
205 doublereal *epstrc,
206 integer *ncfnew);
207static
208int mmtrpj6_(integer *ncofmx,
209 integer *ndimen,
210 integer *ncoeff,
211 doublereal *epsi3d,
212 doublereal *crvlgd,
213 doublereal *ycvmax,
214 doublereal *epstrc,
215 integer *ncfnew);
216static
217integer pow__ii(integer *x,
218 integer *n);
219
220static
221int mvcvin2_(integer *ncoeff,
222 doublereal *crvold,
223 doublereal *crvnew,
224 integer *iercod);
225
226static
227int mvcvinv_(integer *ncoeff,
228 doublereal *crvold,
229 doublereal *crvnew,
230 integer *iercod);
231
232static
233int mvgaus0_(integer *kindic,
234 doublereal *urootl,
235 doublereal *hiltab,
236 integer *nbrval,
237 integer *iercod);
238static
239int mvpscr2_(integer *ncoeff,
240 doublereal *curve2,
241 doublereal *tparam,
242 doublereal *pntcrb);
243
244static
245int mvpscr3_(integer *ncoeff,
246 doublereal *curve2,
247 doublereal *tparam,
248 doublereal *pntcrb);
249
250static struct {
251 doublereal eps1, eps2, eps3, eps4;
252 integer niterm, niterr;
253} mmprcsn_;
254
255static struct {
256 doublereal tdebut, tfinal, verifi, cmherm[576];
257} mmcmher_;
258
259//=======================================================================
260//function : AdvApp2Var_MathBase::mdsptpt_
261//purpose :
262//=======================================================================
263int AdvApp2Var_MathBase::mdsptpt_(integer *ndimen,
264 doublereal *point1,
265 doublereal *point2,
266 doublereal *distan)
267
268{
1ef32e96 269 integer c__8 = 8;
7fd59977 270 /* System generated locals */
271 integer i__1;
272 doublereal d__1;
41194117 273
7fd59977 274 /* Local variables */
1ef32e96
RL
275 integer i__;
276 doublereal* differ = 0;
277 integer ier;
fadcea2c 278 intptr_t iofset, j;
7fd59977 279
280/* **********************************************************************
281*/
282
0d969553 283/* FUNCTION : */
7fd59977 284/* ---------- */
0d969553 285/* CALCULATE DISTANCE BETWEEN TWO POINTS */
7fd59977 286
0d969553 287/* KEYWORDS : */
7fd59977 288/* ----------- */
289/* DISTANCE,POINT. */
290
0d969553 291/* INPUT ARGUMENTS : */
7fd59977 292/* ------------------ */
0d969553
Y
293/* NDIMEN: Space Dimension. */
294/* POINT1: Table of coordinates of the 1st point. */
295/* POINT2: Table of coordinates of the 2nd point. */
7fd59977 296
0d969553 297/* OUTPUT ARGUMENTS : */
7fd59977 298/* ------------------- */
0d969553 299/* DISTAN: Distance between 2 points. */
7fd59977 300
0d969553 301/* COMMONS USED : */
7fd59977 302/* ---------------- */
303
0d969553 304/* REFERENCES CALLED : */
7fd59977 305/* ----------------------- */
306
0d969553 307/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 308/* ----------------------------------- */
7fd59977 309/* > */
310/* **********************************************************************
311*/
312
313
314/* ***********************************************************************
315 */
0d969553 316/* INITIALISATION */
7fd59977 317/* ***********************************************************************
318 */
319
0d969553 320 /* Parameter adjustment */
7fd59977 321 --point2;
322 --point1;
323
324 /* Function Body */
325 iofset = 0;
326 ier = 0;
327
328/* ***********************************************************************
329 */
330/* TRAITEMENT */
331/* ***********************************************************************
332 */
333
1ef32e96 334 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
7fd59977 335 if (*ndimen > 100) {
1ef32e96 336 anAdvApp2Var_SysBase.mcrrqst_(&c__8, ndimen, differ, &iofset, &ier);
7fd59977 337 }
338
0d969553 339/* --- If allocation is refused, the trivial method is applied. */
7fd59977 340
341 if (ier > 0) {
342
343 *distan = 0.;
344 i__1 = *ndimen;
345 for (i__ = 1; i__ <= i__1; ++i__) {
346/* Computing 2nd power */
347 d__1 = point1[i__] - point2[i__];
348 *distan += d__1 * d__1;
349 }
350 *distan = sqrt(*distan);
351
0d969553 352/* --- Otherwise MZSNORM is used to minimize the risks of overflow
7fd59977 353*/
354
355 } else {
356 i__1 = *ndimen;
357 for (i__ = 1; i__ <= i__1; ++i__) {
358 j=iofset + i__ - 1;
359 differ[j] = point2[i__] - point1[i__];
360 }
361
362 *distan = AdvApp2Var_MathBase::mzsnorm_(ndimen, &differ[iofset]);
363
364 }
365
366/* ***********************************************************************
367 */
0d969553 368/* RETURN CALLING PROGRAM */
7fd59977 369/* ***********************************************************************
370 */
371
0d969553 372/* --- Dynamic Desallocation */
7fd59977 373
374 if (iofset != 0) {
1ef32e96 375 anAdvApp2Var_SysBase.mcrdelt_(&c__8, ndimen, differ, &iofset, &ier);
7fd59977 376 }
377
378 return 0 ;
379} /* mdsptpt_ */
380
381//=======================================================================
382//function : mfac_
383//purpose :
384//=======================================================================
385int mfac_(doublereal *f,
386 integer *n)
387
388{
389 /* System generated locals */
390 integer i__1;
391
392 /* Local variables */
1ef32e96 393 integer i__;
7fd59977 394
395/* FORTRAN CONFORME AU TEXT */
396/* CALCUL DE MFACTORIEL N */
397 /* Parameter adjustments */
398 --f;
399
400 /* Function Body */
401 f[1] = (float)1.;
402 i__1 = *n;
403 for (i__ = 2; i__ <= i__1; ++i__) {
404/* L10: */
405 f[i__] = i__ * f[i__ - 1];
406 }
407 return 0;
408} /* mfac_ */
409
410//=======================================================================
411//function : AdvApp2Var_MathBase::mmapcmp_
412//purpose :
413//=======================================================================
414int AdvApp2Var_MathBase::mmapcmp_(integer *ndim,
415 integer *ncofmx,
416 integer *ncoeff,
417 doublereal *crvold,
418 doublereal *crvnew)
419
420{
421 /* System generated locals */
422 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
423 i__2;
41194117 424
7fd59977 425 /* Local variables */
1ef32e96 426 integer ipair, nd, ndegre, impair, ibb, idg;
7fd59977 427 //extern int mgsomsg_();//mgenmsg_(),
7fd59977 428
429/* **********************************************************************
430*/
431
0d969553 432/* FUNCTION : */
7fd59977 433/* ---------- */
0d969553
Y
434/* Compression of curve CRVOLD in a table of */
435/* coeff. of even : CRVNEW(*,0,*) */
436/* and uneven range : CRVNEW(*,1,*). */
7fd59977 437
0d969553 438/* KEYWORDS : */
7fd59977 439/* ----------- */
0d969553 440/* COMPRESSION,CURVE. */
7fd59977 441
0d969553 442/* INPUT ARGUMENTS : */
7fd59977 443/* ------------------ */
0d969553
Y
444/* NDIM : Space Dimension. */
445/* NCOFMX : Max nb of coeff. of the curve to compress. */
446/* NCOEFF : Max nb of coeff. of the compressed curve. */
447/* CRVOLD : The curve (0:NCOFMX-1,NDIM) to compress. */
7fd59977 448
0d969553 449/* OUTPUT ARGUMENTS : */
7fd59977 450/* ------------------- */
0d969553 451/* CRVNEW : Curve compacted in (0:(NCOEFF-1)/2,0,NDIM) (containing
7fd59977 452*/
0d969553
Y
453/* even terms) and in (0:(NCOEFF-1)/2,1,NDIM) */
454/* (containing uneven terms). */
7fd59977 455
0d969553 456/* COMMONS USED : */
7fd59977 457/* ---------------- */
458
0d969553 459/* REFERENCES CALLED : */
7fd59977 460/* ----------------------- */
461
0d969553 462/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 463/* ----------------------------------- */
0d969553
Y
464/* This routine is useful to prepare coefficients of a */
465/* curve in an orthogonal base (Legendre or Jacobi) before */
466/* calculating the coefficients in the canonical; base [-1,1] by */
7fd59977 467/* MMJACAN. */
7fd59977 468/* ***********************************************************************
469 */
470
0d969553 471/* Name of the routine */
7fd59977 472
473 /* Parameter adjustments */
474 crvold_dim1 = *ncofmx;
475 crvold_offset = crvold_dim1;
476 crvold -= crvold_offset;
477 crvnew_dim1 = (*ncoeff - 1) / 2 + 1;
478 crvnew_offset = crvnew_dim1 << 1;
479 crvnew -= crvnew_offset;
480
481 /* Function Body */
482 ibb = AdvApp2Var_SysBase::mnfndeb_();
483 if (ibb >= 3) {
484 AdvApp2Var_SysBase::mgenmsg_("MMAPCMP", 7L);
485 }
486
487 ndegre = *ncoeff - 1;
488 i__1 = *ndim;
489 for (nd = 1; nd <= i__1; ++nd) {
490 ipair = 0;
491 i__2 = ndegre / 2;
492 for (idg = 0; idg <= i__2; ++idg) {
493 crvnew[idg + (nd << 1) * crvnew_dim1] = crvold[ipair + nd *
494 crvold_dim1];
495 ipair += 2;
496/* L200: */
497 }
498 if (ndegre < 1) {
499 goto L400;
500 }
501 impair = 1;
502 i__2 = (ndegre - 1) / 2;
503 for (idg = 0; idg <= i__2; ++idg) {
504 crvnew[idg + ((nd << 1) + 1) * crvnew_dim1] = crvold[impair + nd *
505 crvold_dim1];
506 impair += 2;
507/* L300: */
508 }
509
510L400:
511/* L100: */
512 ;
513 }
514
515/* ---------------------------------- The end ---------------------------
516*/
517
518 if (ibb >= 3) {
519 AdvApp2Var_SysBase::mgsomsg_("MMAPCMP", 7L);
520 }
521 return 0;
522} /* mmapcmp_ */
523
524//=======================================================================
525//function : mmaper0_
526//purpose :
527//=======================================================================
528int mmaper0_(integer *ncofmx,
529 integer *ndimen,
530 integer *ncoeff,
531 doublereal *crvlgd,
532 integer *ncfnew,
533 doublereal *ycvmax,
534 doublereal *errmax)
535
536{
537 /* System generated locals */
538 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
539 doublereal d__1;
41194117 540
7fd59977 541 /* Local variables */
1ef32e96
RL
542 integer ncut;
543 doublereal bidon;
544 integer ii, nd;
7fd59977 545
546/* ***********************************************************************
547 */
548
0d969553 549/* FUNCTION : */
7fd59977 550/* ---------- */
0d969553
Y
551/* Calculate the max error of approximation done when */
552/* only the first NCFNEW coefficients of a curve are preserved.
7fd59977 553*/
0d969553
Y
554/* Degree NCOEFF-1 written in the base of Legendre (Jacobi */
555/* of order 0). */
7fd59977 556
0d969553 557/* KEYWORDS : */
7fd59977 558/* ----------- */
0d969553 559/* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
7fd59977 560
0d969553 561/* INPUT ARGUMENTS : */
7fd59977 562/* ------------------ */
0d969553
Y
563/* NCOFMX : Max. degree of the curve. */
564/* NDIMEN : Space dimension. */
565/* NCOEFF : Degree +1 of the curve. */
566/* CRVLGD : Curve the degree which of should be lowered. */
567/* NCFNEW : Degree +1 of the resulting polynom. */
7fd59977 568
0d969553 569/* OUTPUT ARGUMENTS : */
7fd59977 570/* ------------------- */
0d969553 571/* YCVMAX : Auxiliary Table (max error on each dimension).
7fd59977 572*/
0d969553 573/* ERRMAX : Precision of the approximation. */
7fd59977 574
0d969553 575/* COMMONS USED : */
7fd59977 576/* ---------------- */
577
0d969553 578/* REFERENCES CALLED : */
7fd59977 579/* ----------------------- */
580
0d969553 581/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 582/* ----------------------------------- */
7fd59977 583/* ***********************************************************************
584 */
585
586
0d969553 587/* ------------------- Init to calculate an error -----------------------
7fd59977 588*/
589
590 /* Parameter adjustments */
591 --ycvmax;
592 crvlgd_dim1 = *ncofmx;
593 crvlgd_offset = crvlgd_dim1 + 1;
594 crvlgd -= crvlgd_offset;
595
596 /* Function Body */
597 i__1 = *ndimen;
598 for (ii = 1; ii <= i__1; ++ii) {
599 ycvmax[ii] = 0.;
600/* L100: */
601 }
602
0d969553 603/* ------ Minimum that can be reached : Stop at 1 or NCFNEW ------
7fd59977 604*/
605
606 ncut = 1;
607 if (*ncfnew + 1 > ncut) {
608 ncut = *ncfnew + 1;
609 }
610
0d969553 611/* -------------- Elimination of high degree coefficients-----------
7fd59977 612*/
0d969553 613/* ----------- Loop on the series of Legendre: NCUT --> NCOEFF --------
7fd59977 614*/
615
616 i__1 = *ncoeff;
617 for (ii = ncut; ii <= i__1; ++ii) {
0d969553 618/* Factor of renormalization (Maximum of Li(t)). */
7fd59977 619 bidon = ((ii - 1) * 2. + 1.) / 2.;
620 bidon = sqrt(bidon);
621
622 i__2 = *ndimen;
623 for (nd = 1; nd <= i__2; ++nd) {
41194117 624 ycvmax[nd] += (d__1 = crvlgd[ii + nd * crvlgd_dim1], advapp_abs(d__1)) *
7fd59977 625 bidon;
626/* L310: */
627 }
628/* L300: */
629 }
630
0d969553 631/* -------------- The error is the norm of the vector error ---------------
7fd59977 632*/
633
634 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
635
636/* --------------------------------- Fin --------------------------------
637*/
638
639 return 0;
640} /* mmaper0_ */
641
642//=======================================================================
643//function : mmaper2_
644//purpose :
645//=======================================================================
646int mmaper2_(integer *ncofmx,
647 integer *ndimen,
648 integer *ncoeff,
649 doublereal *crvjac,
650 integer *ncfnew,
651 doublereal *ycvmax,
652 doublereal *errmax)
653
654{
655 /* Initialized data */
656
657 static doublereal xmaxj[57] = { .9682458365518542212948163499456,
658 .986013297183269340427888048593603,
659 1.07810420343739860362585159028115,
660 1.17325804490920057010925920756025,
661 1.26476561266905634732910520370741,
662 1.35169950227289626684434056681946,
663 1.43424378958284137759129885012494,
664 1.51281316274895465689402798226634,
665 1.5878364329591908800533936587012,
666 1.65970112228228167018443636171226,
667 1.72874345388622461848433443013543,
668 1.7952515611463877544077632304216,
669 1.85947199025328260370244491818047,
670 1.92161634324190018916351663207101,
671 1.98186713586472025397859895825157,
672 2.04038269834980146276967984252188,
673 2.09730119173852573441223706382076,
674 2.15274387655763462685970799663412,
675 2.20681777186342079455059961912859,
676 2.25961782459354604684402726624239,
677 2.31122868752403808176824020121524,
678 2.36172618435386566570998793688131,
679 2.41117852396114589446497298177554,
680 2.45964731268663657873849811095449,
681 2.50718840313973523778244737914028,
682 2.55385260994795361951813645784034,
683 2.59968631659221867834697883938297,
684 2.64473199258285846332860663371298,
685 2.68902863641518586789566216064557,
686 2.73261215675199397407027673053895,
687 2.77551570192374483822124304745691,
688 2.8177699459714315371037628127545,
689 2.85940333797200948896046563785957,
690 2.90044232019793636101516293333324,
691 2.94091151970640874812265419871976,
692 2.98083391718088702956696303389061,
693 3.02023099621926980436221568258656,
694 3.05912287574998661724731962377847,
695 3.09752842783622025614245706196447,
696 3.13546538278134559341444834866301,
697 3.17295042316122606504398054547289,
698 3.2099992681699613513775259670214,
699 3.24662674946606137764916854570219,
700 3.28284687953866689817670991319787,
701 3.31867291347259485044591136879087,
702 3.35411740487202127264475726990106,
703 3.38919225660177218727305224515862,
704 3.42390876691942143189170489271753,
705 3.45827767149820230182596660024454,
706 3.49230918177808483937957161007792,
707 3.5260130200285724149540352829756,
708 3.55939845146044235497103883695448,
709 3.59247431368364585025958062194665,
710 3.62524904377393592090180712976368,
711 3.65773070318071087226169680450936,
712 3.68992700068237648299565823810245,
713 3.72184531357268220291630708234186 };
714
715 /* System generated locals */
716 integer crvjac_dim1, crvjac_offset, i__1, i__2;
717 doublereal d__1;
718
719 /* Local variables */
1ef32e96
RL
720 integer idec, ncut;
721 doublereal bidon;
722 integer ii, nd;
7fd59977 723
724
725
726/* ***********************************************************************
727 */
728
729/* FONCTION : */
730/* ---------- */
0d969553 731/* Calculate max approximation error i faite lorsque l' on */
7fd59977 732/* ne conserve que les premiers NCFNEW coefficients d' une courbe
733*/
734/* de degre NCOEFF-1 ecrite dans la base de Jacobi d' ordre 2. */
735
0d969553 736/* KEYWORDS : */
7fd59977 737/* ----------- */
0d969553
Y
738/* JACOBI, POLYGON, APPROXIMATION, ERROR. */
739/*
740/* INPUT ARGUMENTS : */
7fd59977 741/* ------------------ */
0d969553
Y
742/* NCOFMX : Max. degree of the curve. */
743/* NDIMEN : Space dimension. */
744/* NCOEFF : Degree +1 of the curve. */
745/* CRVLGD : Curve the degree which of should be lowered. */
746/* NCFNEW : Degree +1 of the resulting polynom. */
7fd59977 747
0d969553 748/* OUTPUT ARGUMENTS : */
7fd59977 749/* ------------------- */
0d969553 750/* YCVMAX : Auxiliary Table (max error on each dimension).
7fd59977 751*/
0d969553 752/* ERRMAX : Precision of the approximation. */
7fd59977 753
0d969553 754/* COMMONS USED : */
7fd59977 755/* ---------------- */
756
0d969553 757/* REFERENCES CALLED : */
7fd59977 758/* ----------------------- */
0d969553 759/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 760/* ----------------------------------- */
761
7fd59977 762
763
0d969553 764/* ------------------ Table of maximums of (1-t2)*Ji(t) ----------------
7fd59977 765*/
766
767 /* Parameter adjustments */
768 --ycvmax;
769 crvjac_dim1 = *ncofmx;
770 crvjac_offset = crvjac_dim1 + 1;
771 crvjac -= crvjac_offset;
772
773 /* Function Body */
774
775
776
0d969553 777/* ------------------- Init for error calculation -----------------------
7fd59977 778*/
779
780 i__1 = *ndimen;
781 for (ii = 1; ii <= i__1; ++ii) {
782 ycvmax[ii] = 0.;
783/* L100: */
784 }
785
0d969553 786/* ------ Min. Degree that can be attained : Stop at 3 or NCFNEW ------
7fd59977 787*/
788
789 idec = 3;
790/* Computing MAX */
791 i__1 = idec, i__2 = *ncfnew + 1;
41194117 792 ncut = advapp_max(i__1,i__2);
7fd59977 793
0d969553 794/* -------------- Removal of coefficients of high degree -----------
7fd59977 795*/
0d969553 796/* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
7fd59977 797*/
798
799 i__1 = *ncoeff;
800 for (ii = ncut; ii <= i__1; ++ii) {
0d969553 801/* Factor of renormalization. */
7fd59977 802 bidon = xmaxj[ii - idec];
803 i__2 = *ndimen;
804 for (nd = 1; nd <= i__2; ++nd) {
41194117 805 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
7fd59977 806 bidon;
807/* L310: */
808 }
809/* L300: */
810 }
811
0d969553 812/* -------------- The error is the norm of the vector error ---------------
7fd59977 813*/
814
815 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
816
817/* --------------------------------- Fin --------------------------------
818*/
819
820 return 0;
821} /* mmaper2_ */
822
823/* MAPER4.f -- translated by f2c (version 19960827).
824 You must link the resulting object file with the libraries:
825 -lf2c -lm (in that order)
826*/
827
828/* Subroutine */
829//=======================================================================
830//function : mmaper4_
831//purpose :
832//=======================================================================
833int mmaper4_(integer *ncofmx,
834 integer *ndimen,
835 integer *ncoeff,
836 doublereal *crvjac,
837 integer *ncfnew,
838 doublereal *ycvmax,
839 doublereal *errmax)
840{
841 /* Initialized data */
842
843 static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
844 1.05299572648705464724876659688996,
845 1.0949715351434178709281698645813,
846 1.15078388379719068145021100764647,
847 1.2094863084718701596278219811869,
848 1.26806623151369531323304177532868,
849 1.32549784426476978866302826176202,
850 1.38142537365039019558329304432581,
851 1.43575531950773585146867625840552,
852 1.48850442653629641402403231015299,
853 1.53973611681876234549146350844736,
854 1.58953193485272191557448229046492,
855 1.63797820416306624705258190017418,
856 1.68515974143594899185621942934906,
857 1.73115699602477936547107755854868,
858 1.77604489805513552087086912113251,
859 1.81989256661534438347398400420601,
860 1.86276344480103110090865609776681,
861 1.90471563564740808542244678597105,
862 1.94580231994751044968731427898046,
863 1.98607219357764450634552790950067,
864 2.02556989246317857340333585562678,
865 2.06433638992049685189059517340452,
866 2.10240936014742726236706004607473,
867 2.13982350649113222745523925190532,
868 2.17661085564771614285379929798896,
869 2.21280102016879766322589373557048,
870 2.2484214321456956597803794333791,
871 2.28349755104077956674135810027654,
872 2.31805304852593774867640120860446,
873 2.35210997297725685169643559615022,
874 2.38568889602346315560143377261814,
875 2.41880904328694215730192284109322,
876 2.45148841120796359750021227795539,
877 2.48374387161372199992570528025315,
878 2.5155912654873773953959098501893,
879 2.54704548720896557684101746505398,
880 2.57812056037881628390134077704127,
881 2.60882970619319538196517982945269,
882 2.63918540521920497868347679257107,
883 2.66919945330942891495458446613851,
884 2.69888301230439621709803756505788,
885 2.72824665609081486737132853370048,
886 2.75730041251405791603760003778285,
887 2.78605380158311346185098508516203,
888 2.81451587035387403267676338931454,
889 2.84269522483114290814009184272637,
890 2.87060005919012917988363332454033,
891 2.89823818258367657739520912946934,
892 2.92561704377132528239806135133273,
893 2.95274375377994262301217318010209,
894 2.97962510678256471794289060402033,
895 3.00626759936182712291041810228171,
896 3.03267744830655121818899164295959,
897 3.05886060707437081434964933864149 };
898
899 /* System generated locals */
900 integer crvjac_dim1, crvjac_offset, i__1, i__2;
901 doublereal d__1;
902
903 /* Local variables */
1ef32e96
RL
904 integer idec, ncut;
905 doublereal bidon;
906 integer ii, nd;
7fd59977 907
908
909
910/* ***********************************************************************
911 */
912
0d969553 913/* FUNCTION : */
7fd59977 914/* ---------- */
0d969553
Y
915/* Calculate the max. error of approximation made when */
916/* only first NCFNEW coefficients of a curve are preserved
7fd59977 917*/
0d969553
Y
918/* degree NCOEFF-1 is written in the base of Jacobi of order 4. */
919/* KEYWORDS : */
7fd59977 920/* ----------- */
0d969553 921/* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
7fd59977 922
0d969553 923/* INPUT ARGUMENTS : */
7fd59977 924/* ------------------ */
0d969553
Y
925/* NCOFMX : Max. degree of the curve. */
926/* NDIMEN : Space dimension. */
927/* NCOEFF : Degree +1 of the curve. */
928/* CRVJAC : Curve the degree which of should be lowered. */
929/* NCFNEW : Degree +1 of the resulting polynom. */
7fd59977 930
0d969553 931/* OUTPUT ARGUMENTS : */
7fd59977 932/* ------------------- */
0d969553 933/* YCVMAX : Auxiliary Table (max error on each dimension).
7fd59977 934*/
0d969553 935/* ERRMAX : Precision of the approximation. */
7fd59977 936
0d969553 937/* COMMONS USED : */
7fd59977 938/* ---------------- */
939
0d969553 940/* REFERENCES CALLED : */
7fd59977 941/* ----------------------- */
942
0d969553
Y
943/* DESCRIPTION/NOTES/LIMITATIONS : */
944
7fd59977 945
7fd59977 946/* ***********************************************************************
947 */
948
949
0d969553 950/* ---------------- Table of maximums of ((1-t2)2)*Ji(t) ---------------
7fd59977 951*/
952
953 /* Parameter adjustments */
954 --ycvmax;
955 crvjac_dim1 = *ncofmx;
956 crvjac_offset = crvjac_dim1 + 1;
957 crvjac -= crvjac_offset;
958
959 /* Function Body */
960
961
962
0d969553 963/* ------------------- Init for error calculation -----------------------
7fd59977 964*/
965
966 i__1 = *ndimen;
967 for (ii = 1; ii <= i__1; ++ii) {
968 ycvmax[ii] = 0.;
969/* L100: */
970 }
971
0d969553 972/* ------ Min. Degree that can be attained : Stop at 5 or NCFNEW ------
7fd59977 973*/
974
975 idec = 5;
976/* Computing MAX */
977 i__1 = idec, i__2 = *ncfnew + 1;
41194117 978 ncut = advapp_max(i__1,i__2);
7fd59977 979
0d969553 980/* -------------- Removal of high degree coefficients -----------
7fd59977 981*/
0d969553 982/* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
7fd59977 983*/
984
985 i__1 = *ncoeff;
986 for (ii = ncut; ii <= i__1; ++ii) {
0d969553 987/* Factor of renormalisation. */
7fd59977 988 bidon = xmaxj[ii - idec];
989 i__2 = *ndimen;
990 for (nd = 1; nd <= i__2; ++nd) {
41194117 991 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
7fd59977 992 bidon;
993/* L310: */
994 }
995/* L300: */
996 }
997
0d969553 998/* -------------- The error is the norm of the error vector ---------------
7fd59977 999*/
1000
1001 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
1002
0d969553 1003/* --------------------------------- End --------------------------------
7fd59977 1004*/
1005
1006 return 0;
1007} /* mmaper4_ */
1008
1009//=======================================================================
1010//function : mmaper6_
1011//purpose :
1012//=======================================================================
1013int mmaper6_(integer *ncofmx,
1014 integer *ndimen,
1015 integer *ncoeff,
1016 doublereal *crvjac,
1017 integer *ncfnew,
1018 doublereal *ycvmax,
1019 doublereal *errmax)
1020
1021{
1022 /* Initialized data */
1023
1024 static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
1025 1.11626917091567929907256116528817,
1026 1.1327140810290884106278510474203,
1027 1.1679452722668028753522098022171,
1028 1.20910611986279066645602153641334,
1029 1.25228283758701572089625983127043,
1030 1.29591971597287895911380446311508,
1031 1.3393138157481884258308028584917,
1032 1.3821288728999671920677617491385,
1033 1.42420414683357356104823573391816,
1034 1.46546895108549501306970087318319,
1035 1.50590085198398789708599726315869,
1036 1.54550385142820987194251585145013,
1037 1.58429644271680300005206185490937,
1038 1.62230484071440103826322971668038,
1039 1.65955905239130512405565733793667,
1040 1.69609056468292429853775667485212,
1041 1.73193098017228915881592458573809,
1042 1.7671112206990325429863426635397,
1043 1.80166107681586964987277458875667,
1044 1.83560897003644959204940535551721,
1045 1.86898184653271388435058371983316,
1046 1.90180515174518670797686768515502,
1047 1.93410285411785808749237200054739,
1048 1.96589749778987993293150856865539,
1049 1.99721027139062501070081653790635,
1050 2.02806108474738744005306947877164,
1051 2.05846864831762572089033752595401,
1052 2.08845055210580131460156962214748,
1053 2.11802334209486194329576724042253,
1054 2.14720259305166593214642386780469,
1055 2.17600297710595096918495785742803,
1056 2.20443832785205516555772788192013,
1057 2.2325216999457379530416998244706,
1058 2.2602654243075083168599953074345,
1059 2.28768115912702794202525264301585,
1060 2.3147799369092684021274946755348,
1061 2.34157220782483457076721300512406,
1062 2.36806787963276257263034969490066,
1063 2.39427635443992520016789041085844,
1064 2.42020656255081863955040620243062,
1065 2.44586699364757383088888037359254,
1066 2.47126572552427660024678584642791,
1067 2.49641045058324178349347438430311,
1068 2.52130850028451113942299097584818,
1069 2.54596686772399937214920135190177,
1070 2.5703922285006754089328998222275,
1071 2.59459096001908861492582631591134,
1072 2.61856915936049852435394597597773,
1073 2.64233265984385295286445444361827,
1074 2.66588704638685848486056711408168,
1075 2.68923766976735295746679957665724,
1076 2.71238965987606292679677228666411 };
1077
1078 /* System generated locals */
1079 integer crvjac_dim1, crvjac_offset, i__1, i__2;
1080 doublereal d__1;
1081
1082 /* Local variables */
1ef32e96
RL
1083 integer idec, ncut;
1084 doublereal bidon;
1085 integer ii, nd;
7fd59977 1086
1087
1088
1089/* ***********************************************************************
1090 */
0d969553 1091/* FUNCTION : */
7fd59977 1092/* ---------- */
0d969553
Y
1093/* Calculate the max. error of approximation made when */
1094/* only first NCFNEW coefficients of a curve are preserved
7fd59977 1095*/
0d969553
Y
1096/* degree NCOEFF-1 is written in the base of Jacobi of order 6. */
1097/* KEYWORDS : */
7fd59977 1098/* ----------- */
0d969553 1099/* JACOBI,POLYGON,APPROXIMATION,ERROR. */
7fd59977 1100
0d969553 1101/* INPUT ARGUMENTS : */
7fd59977 1102/* ------------------ */
0d969553
Y
1103/* NCOFMX : Max. degree of the curve. */
1104/* NDIMEN : Space dimension. */
1105/* NCOEFF : Degree +1 of the curve. */
1106/* CRVJAC : Curve the degree which of should be lowered. */
1107/* NCFNEW : Degree +1 of the resulting polynom. */
7fd59977 1108
0d969553 1109/* OUTPUT ARGUMENTS : */
7fd59977 1110/* ------------------- */
0d969553 1111/* YCVMAX : Auxiliary Table (max error on each dimension).
7fd59977 1112*/
0d969553 1113/* ERRMAX : Precision of the approximation. */
7fd59977 1114
0d969553 1115/* COMMONS USED : */
7fd59977 1116/* ---------------- */
1117
0d969553 1118/* REFERENCES CALLED : */
7fd59977 1119/* ----------------------- */
1120
0d969553 1121/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 1122/* > */
1123/* ***********************************************************************
1124 */
1125
1126
0d969553 1127/* ---------------- Table of maximums of ((1-t2)3)*Ji(t) ---------------
7fd59977 1128*/
1129
1130 /* Parameter adjustments */
1131 --ycvmax;
1132 crvjac_dim1 = *ncofmx;
1133 crvjac_offset = crvjac_dim1 + 1;
1134 crvjac -= crvjac_offset;
1135
1136 /* Function Body */
1137
1138
1139
0d969553 1140/* ------------------- Init for error calculation -----------------------
7fd59977 1141*/
1142
1143 i__1 = *ndimen;
1144 for (ii = 1; ii <= i__1; ++ii) {
1145 ycvmax[ii] = 0.;
1146/* L100: */
1147 }
1148
0d969553 1149/* ------ Min Degree that can be attained : Stop at 3 or NCFNEW ------
7fd59977 1150*/
1151
1152 idec = 7;
1153/* Computing MAX */
1154 i__1 = idec, i__2 = *ncfnew + 1;
41194117 1155 ncut = advapp_max(i__1,i__2);
7fd59977 1156
0d969553 1157/* -------------- Removal of high degree coefficients -----------
7fd59977 1158*/
0d969553 1159/* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
7fd59977 1160*/
1161
1162 i__1 = *ncoeff;
1163 for (ii = ncut; ii <= i__1; ++ii) {
0d969553 1164/* Factor of renormalization. */
7fd59977 1165 bidon = xmaxj[ii - idec];
1166 i__2 = *ndimen;
1167 for (nd = 1; nd <= i__2; ++nd) {
41194117 1168 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
7fd59977 1169 bidon;
1170/* L310: */
1171 }
1172/* L300: */
1173 }
1174
0d969553 1175/* -------------- The error is the norm of the vector error ---------------
7fd59977 1176*/
1177
1178 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
1179
0d969553 1180/* --------------------------------- END --------------------------------
7fd59977 1181*/
1182
1183 return 0;
1184} /* mmaper6_ */
1185
1186//=======================================================================
1187//function : AdvApp2Var_MathBase::mmaperx_
1188//purpose :
1189//=======================================================================
1190int AdvApp2Var_MathBase::mmaperx_(integer *ncofmx,
1191 integer *ndimen,
1192 integer *ncoeff,
1193 integer *iordre,
1194 doublereal *crvjac,
1195 integer *ncfnew,
1196 doublereal *ycvmax,
1197 doublereal *errmax,
1198 integer *iercod)
1199
1200{
1201 /* System generated locals */
1202 integer crvjac_dim1, crvjac_offset;
41194117 1203
7fd59977 1204 /* Local variables */
1ef32e96 1205 integer jord;
7fd59977 1206
1207/* **********************************************************************
1208*/
0d969553 1209/* FUNCTION : */
7fd59977 1210/* ---------- */
0d969553
Y
1211/* Calculate the max. error of approximation made when */
1212/* only first NCFNEW coefficients of a curve are preserved
7fd59977 1213*/
0d969553
Y
1214/* degree NCOEFF-1 is written in the base of Jacobi of order IORDRE. */
1215/* KEYWORDS : */
7fd59977 1216/* ----------- */
0d969553 1217/* JACOBI,LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
7fd59977 1218
0d969553 1219/* INPUT ARGUMENTS : */
7fd59977 1220/* ------------------ */
0d969553
Y
1221/* NCOFMX : Max. degree of the curve. */
1222/* NDIMEN : Space dimension. */
1223/* NCOEFF : Degree +1 of the curve. */
1224/* IORDRE : Order of continuity at the extremities. */
1225/* CRVJAC : Curve the degree which of should be lowered. */
1226/* NCFNEW : Degree +1 of the resulting polynom. */
1227
1228/* OUTPUT ARGUMENTS : */
7fd59977 1229/* ------------------- */
0d969553
Y
1230/* YCVMAX : Auxiliary Table (max error on each dimension).
1231*/
1232/* ERRMAX : Precision of the approximation. */
7fd59977 1233/* IERCOD = 0, OK */
0d969553
Y
1234/* = 1, order of constraints (IORDRE) is not within the */
1235/* autorized values. */
1236/* COMMONS USED : */
7fd59977 1237/* ---------------- */
1238
0d969553 1239/* REFERENCES CALLED : */
7fd59977 1240/* ----------------------- */
1241
0d969553 1242/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 1243/* ----------------------------------- */
0d969553 1244/* Canceled and replaced MMAPERR. */
7fd59977 1245/* ***********************************************************************
1246 */
1247
1248
1249 /* Parameter adjustments */
1250 --ycvmax;
1251 crvjac_dim1 = *ncofmx;
1252 crvjac_offset = crvjac_dim1 + 1;
1253 crvjac -= crvjac_offset;
1254
1255 /* Function Body */
1256 *iercod = 0;
0d969553 1257/* --> Order of Jacobi polynoms */
7fd59977 1258 jord = ( *iordre + 1) << 1;
1259
1260 if (jord == 0) {
1261 mmaper0_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1262 ycvmax[1], errmax);
1263 } else if (jord == 2) {
1264 mmaper2_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1265 ycvmax[1], errmax);
1266 } else if (jord == 4) {
1267 mmaper4_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1268 ycvmax[1], errmax);
1269 } else if (jord == 6) {
1270 mmaper6_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1271 ycvmax[1], errmax);
1272 } else {
1273 *iercod = 1;
1274 }
1275
1276/* ----------------------------------- Fin ------------------------------
1277*/
1278
1279 return 0;
1280} /* mmaperx_ */
1281
1282//=======================================================================
1283//function : mmarc41_
1284//purpose :
1285//=======================================================================
1286 int mmarc41_(integer *ndimax,
1287 integer *ndimen,
1288 integer *ncoeff,
1289 doublereal *crvold,
1290 doublereal *upara0,
1291 doublereal *upara1,
1292 doublereal *crvnew,
1293 integer *iercod)
1294
1295{
1296 /* System generated locals */
1297 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
1298 i__2, i__3;
1299
1300 /* Local variables */
1ef32e96
RL
1301 integer nboct;
1302 doublereal tbaux[61];
1303 integer nd;
1304 doublereal bid;
1305 integer ncf, ncj;
7fd59977 1306
1307
1308/* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
1309/* IMPLICIT INTEGER (I-N) */
1310
1311/* ***********************************************************************
1312 */
1313
0d969553 1314/* FUNCTION : */
7fd59977 1315/* ---------- */
0d969553
Y
1316/* Creation of curve C2(v) defined on (0,1) identic to */
1317/* curve C1(u) defined on (U0,U1) (change of parameter */
1318/* of a curve). */
7fd59977 1319
0d969553 1320/* KEYWORDS : */
7fd59977 1321/* ----------- */
0d969553 1322/* LIMITATION, RESTRICTION, CURVE */
7fd59977 1323
0d969553 1324/* INPUT ARGUMENTS : */
7fd59977 1325/* ------------------ */
0d969553
Y
1326/* NDIMAX : Space Dimensioning. */
1327/* NDIMEN : Curve Dimension. */
1328/* NCOEFF : Nb of coefficients of the curve. */
1329/* CRVOLD : Curve to be limited. */
1330/* UPARA0 : Min limit of the interval limiting the curve.
7fd59977 1331*/
0d969553 1332/* UPARA1 : Max limit of the interval limiting the curve.
7fd59977 1333*/
1334
0d969553 1335/* OUTPUT ARGUMENTS : */
7fd59977 1336/* ------------------- */
0d969553
Y
1337/* CRVNEW : Relimited curve, defined on (0,1) and equal to */
1338/* CRVOLD defined on (U0,U1). */
7fd59977 1339/* IERCOD : = 0, OK */
0d969553 1340/* =10, Nb of coeff. <1 or > 61. */
7fd59977 1341
0d969553 1342/* COMMONS USED : */
7fd59977 1343/* ---------------- */
0d969553 1344/* REFERENCES CALLED : */
7fd59977 1345/* ---------------------- */
1346/* Type Name */
1347/* MAERMSG MCRFILL MVCVIN2 */
1348/* MVCVINV */
1349
0d969553 1350/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 1351/* ----------------------------------- */
0d969553
Y
1352/* ---> Algorithm used in this general case is based on the */
1353/* following principle : */
1354/* Let S(t) = a0 + a1*t + a2*t**2 + ... of degree NCOEFF-1, and */
1355/* U(t) = b0 + b1*t, then the coeff. of */
1356/* S(U(t)) are calculated step by step with help of table TBAUX. */
1357/* At each step number N (N=2 to NCOEFF), TBAUX(n) contains */
1358/* the n-th coefficient of U(t)**N for n=1 to N. (RBD) */
7fd59977 1359/* ---> Reference : KNUTH, 'The Art of Computer Programming', */
1360/* Vol. 2/'Seminumerical Algorithms', */
1361/* Ex. 11 p:451 et solution p:562. (RBD) */
1362
0d969553
Y
1363/* ---> Removal of the input argument CRVOLD by CRVNEW is */
1364/* possible, which means that the call : */
7fd59977 1365/* CALL MMARC41(NDIMAX,NDIMEN,NCOEFF,CURVE,UPARA0,UPARA1 */
1366/* ,CURVE,IERCOD) */
0d969553 1367/* is absolutely LEGAL. (RBD) */
7fd59977 1368
1369/* > */
1370/* **********************************************************************
1371*/
1372
0d969553 1373/* Name of the routine */
7fd59977 1374
0d969553
Y
1375/* Auxiliary table of coefficients of (UPARA1-UPARA0)T+UPARA0 */
1376/* with power N=1 to NCOEFF-1. */
7fd59977 1377
1378
1379 /* Parameter adjustments */
1380 crvnew_dim1 = *ndimax;
1381 crvnew_offset = crvnew_dim1 + 1;
1382 crvnew -= crvnew_offset;
1383 crvold_dim1 = *ndimax;
1384 crvold_offset = crvold_dim1 + 1;
1385 crvold -= crvold_offset;
1386
1387 /* Function Body */
1388 *iercod = 0;
1389/* **********************************************************************
1390*/
0d969553 1391/* CASE WHEN PROCESSING CAN'T BE DONE */
7fd59977 1392/* **********************************************************************
1393*/
1394 if (*ncoeff > 61 || *ncoeff < 1) {
1395 *iercod = 10;
1396 goto L9999;
1397 }
1398/* **********************************************************************
1399*/
0d969553 1400/* IF NO CHANGES */
7fd59977 1401/* **********************************************************************
1402*/
1403 if (*ndimen == *ndimax && *upara0 == 0. && *upara1 == 1.) {
1404 nboct = (*ndimax << 3) * *ncoeff;
fadcea2c
RL
1405 AdvApp2Var_SysBase::mcrfill_(&nboct,
1406 &crvold[crvold_offset],
1407 &crvnew[crvnew_offset]);
7fd59977 1408 goto L9999;
1409 }
1410/* **********************************************************************
1411*/
0d969553 1412/* INVERSION 3D : FAST PROCESSING */
7fd59977 1413/* **********************************************************************
1414*/
1415 if (*upara0 == 1. && *upara1 == 0.) {
1416 if (*ndimen == 3 && *ndimax == 3 && *ncoeff <= 21) {
1417 mvcvinv_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
1418 iercod);
1419 goto L9999;
1420 }
1421/* ******************************************************************
1422**** */
0d969553 1423/* INVERSION 2D : FAST PROCESSING */
7fd59977 1424/* ******************************************************************
1425**** */
1426 if (*ndimen == 2 && *ndimax == 2 && *ncoeff <= 21) {
1427 mvcvin2_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
1428 iercod);
1429 goto L9999;
1430 }
1431 }
1432/* **********************************************************************
1433*/
0d969553 1434/* GENERAL PROCESSING */
7fd59977 1435/* **********************************************************************
1436*/
0d969553 1437/* -------------------------- Initializations ---------------------------
7fd59977 1438*/
1439
1440 i__1 = *ndimen;
1441 for (nd = 1; nd <= i__1; ++nd) {
1442 crvnew[nd + crvnew_dim1] = crvold[nd + crvold_dim1];
1443/* L100: */
1444 }
1445 if (*ncoeff == 1) {
1446 goto L9999;
1447 }
1448 tbaux[0] = *upara0;
1449 tbaux[1] = *upara1 - *upara0;
1450
0d969553 1451/* ----------------------- Calculation of coeff. of CRVNEW ------------------
7fd59977 1452*/
1453
1454 i__1 = *ncoeff - 1;
1455 for (ncf = 2; ncf <= i__1; ++ncf) {
1456
0d969553 1457/* ------------ Take into account NCF-th coeff. of CRVOLD --------
7fd59977 1458---- */
1459
1460 i__2 = ncf - 1;
1461 for (ncj = 1; ncj <= i__2; ++ncj) {
1462 bid = tbaux[ncj - 1];
1463 i__3 = *ndimen;
1464 for (nd = 1; nd <= i__3; ++nd) {
1465 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + ncf *
1466 crvold_dim1] * bid;
1467/* L400: */
1468 }
1469/* L300: */
1470 }
1471
1472 bid = tbaux[ncf - 1];
1473 i__2 = *ndimen;
1474 for (nd = 1; nd <= i__2; ++nd) {
1475 crvnew[nd + ncf * crvnew_dim1] = crvold[nd + ncf * crvold_dim1] *
1476 bid;
1477/* L500: */
1478 }
1479
0d969553 1480/* --------- Calculate (NCF+1) coeff. of ((U1-U0)*t + U0)**(NCF) ---
7fd59977 1481---- */
1482
1483 bid = *upara1 - *upara0;
1484 tbaux[ncf] = tbaux[ncf - 1] * bid;
1485 for (ncj = ncf; ncj >= 2; --ncj) {
1486 tbaux[ncj - 1] = tbaux[ncj - 1] * *upara0 + tbaux[ncj - 2] * bid;
1487/* L600: */
1488 }
1489 tbaux[0] *= *upara0;
1490
1491/* L200: */
1492 }
1493
0d969553 1494/* -------------- Take into account the last coeff. of CRVOLD -----------
7fd59977 1495*/
1496
1497 i__1 = *ncoeff - 1;
1498 for (ncj = 1; ncj <= i__1; ++ncj) {
1499 bid = tbaux[ncj - 1];
1500 i__2 = *ndimen;
1501 for (nd = 1; nd <= i__2; ++nd) {
1502 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + *ncoeff *
1503 crvold_dim1] * bid;
1504/* L800: */
1505 }
1506/* L700: */
1507 }
1508 i__1 = *ndimen;
1509 for (nd = 1; nd <= i__1; ++nd) {
1510 crvnew[nd + *ncoeff * crvnew_dim1] = crvold[nd + *ncoeff *
1511 crvold_dim1] * tbaux[*ncoeff - 1];
1512/* L900: */
1513 }
1514
1515/* ---------------------------- The end ---------------------------------
1516*/
1517
1518L9999:
1519 if (*iercod != 0) {
1520 AdvApp2Var_SysBase::maermsg_("MMARC41", iercod, 7L);
1521 }
1522
1523 return 0 ;
1524} /* mmarc41_ */
1525
1526//=======================================================================
1527//function : AdvApp2Var_MathBase::mmarcin_
1528//purpose :
1529//=======================================================================
1530int AdvApp2Var_MathBase::mmarcin_(integer *ndimax,
1531 integer *ndim,
1532 integer *ncoeff,
1533 doublereal *crvold,
1534 doublereal *u0,
1535 doublereal *u1,
1536 doublereal *crvnew,
1537 integer *iercod)
1538
1539{
1540 /* System generated locals */
1541 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
1542 i__2, i__3;
1543 doublereal d__1;
1544
1545 /* Local variables */
1ef32e96
RL
1546 doublereal x0, x1;
1547 integer nd;
1548 doublereal tabaux[61];
1549 integer ibb;
1550 doublereal bid;
1551 integer ncf;
1552 integer ncj;
1553 doublereal eps3;
7fd59977 1554
1555
1556
1557/* **********************************************************************
0d969553 1558*//* FUNCTION : */
7fd59977 1559/* ---------- */
0d969553
Y
1560/* Creation of curve C2(v) defined on [U0,U1] identic to */
1561/* curve C1(u) defined on [-1,1] (change of parameter */
1562/* of a curve) with INVERSION of indices of the resulting table. */
7fd59977 1563
0d969553 1564/* KEYWORDS : */
7fd59977 1565/* ----------- */
0d969553 1566/* GENERALIZED LIMITATION, RESTRICTION, INVERSION, CURVE */
7fd59977 1567
0d969553 1568/* INPUT ARGUMENTS : */
7fd59977 1569/* ------------------ */
0d969553
Y
1570/* NDIMAX : Maximum Space Dimensioning. */
1571/* NDIMEN : Curve Dimension. */
1572/* NCOEFF : Nb of coefficients of the curve. */
1573/* CRVOLD : Curve to be limited. */
1574/* U0 : Min limit of the interval limiting the curve.
1575*/
1576/* U1 : Max limit of the interval limiting the curve.
1577*/
1578
1579/* OUTPUT ARGUMENTS : */
7fd59977 1580/* ------------------- */
0d969553
Y
1581/* CRVNEW : Relimited curve, defined on [U0,U1] and equal to */
1582/* CRVOLD defined on [-1,1]. */
7fd59977 1583/* IERCOD : = 0, OK */
0d969553
Y
1584/* =10, Nb of coeff. <1 or > 61. */
1585/* =13, the requested interval of variation is null. */
1586/* COMMONS USED : */
7fd59977 1587/* ---------------- */
0d969553
Y
1588/* REFERENCES CALLED : */
1589/* ---------------------- */
1590/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 1591/* ----------------------------------- */
7fd59977 1592/* > */
1593/* **********************************************************************
1594*/
1595
0d969553 1596/* Name of the routine */
7fd59977 1597
0d969553
Y
1598/* Auxiliary table of coefficients of X1*T+X0 */
1599/* with power N=1 to NCOEFF-1. */
7fd59977 1600
1601
1602 /* Parameter adjustments */
1603 crvnew_dim1 = *ndimax;
1604 crvnew_offset = crvnew_dim1 + 1;
1605 crvnew -= crvnew_offset;
1606 crvold_dim1 = *ncoeff;
1607 crvold_offset = crvold_dim1 + 1;
1608 crvold -= crvold_offset;
1609
1610 /* Function Body */
1611 ibb = AdvApp2Var_SysBase::mnfndeb_();
1612 if (ibb >= 2) {
1613 AdvApp2Var_SysBase::mgenmsg_("MMARCIN", 7L);
1614 }
1615
0d969553 1616/* At zero machine it is tested if the output interval is not null */
7fd59977 1617
1618 AdvApp2Var_MathBase::mmveps3_(&eps3);
41194117 1619 if ((d__1 = *u1 - *u0, advapp_abs(d__1)) < eps3) {
7fd59977 1620 *iercod = 13;
1621 goto L9999;
1622 }
1623 *iercod = 0;
1624
1625/* **********************************************************************
1626*/
0d969553 1627/* CASE WHEN THE PROCESSING IS IMPOSSIBLE */
7fd59977 1628/* **********************************************************************
1629*/
1630 if (*ncoeff > 61 || *ncoeff < 1) {
1631 *iercod = 10;
1632 goto L9999;
1633 }
1634/* **********************************************************************
1635*/
0d969553
Y
1636/* IF NO CHANGE OF THE INTERVAL OF DEFINITION */
1637/* (ONLY INVERSION OF INDICES OF TABLE CRVOLD) */
7fd59977 1638/* **********************************************************************
1639*/
1640 if (*ndim == *ndimax && *u0 == -1. && *u1 == 1.) {
1641 AdvApp2Var_MathBase::mmcvinv_(ndim, ncoeff, ndim, &crvold[crvold_offset], &crvnew[
1642 crvnew_offset]);
1643 goto L9999;
1644 }
1645/* **********************************************************************
1646*/
0d969553 1647/* CASE WHEN THE NEW INTERVAL OF DEFINITION IS [0,1] */
7fd59977 1648/* **********************************************************************
1649*/
1650 if (*u0 == 0. && *u1 == 1.) {
1651 mmcvstd_(ncoeff, ndimax, ncoeff, ndim, &crvold[crvold_offset], &
1652 crvnew[crvnew_offset]);
1653 goto L9999;
1654 }
1655/* **********************************************************************
1656*/
0d969553 1657/* GENERAL PROCESSING */
7fd59977 1658/* **********************************************************************
1659*/
0d969553 1660/* -------------------------- Initialization ---------------------------
7fd59977 1661*/
1662
1663 x0 = -(*u1 + *u0) / (*u1 - *u0);
1664 x1 = 2. / (*u1 - *u0);
1665 i__1 = *ndim;
1666 for (nd = 1; nd <= i__1; ++nd) {
1667 crvnew[nd + crvnew_dim1] = crvold[nd * crvold_dim1 + 1];
1668/* L100: */
1669 }
1670 if (*ncoeff == 1) {
1671 goto L9999;
1672 }
1673 tabaux[0] = x0;
1674 tabaux[1] = x1;
1675
0d969553 1676/* ----------------------- Calculation of coeff. of CRVNEW ------------------
7fd59977 1677*/
1678
1679 i__1 = *ncoeff - 1;
1680 for (ncf = 2; ncf <= i__1; ++ncf) {
1681
0d969553 1682/* ------------ Take into account the NCF-th coeff. of CRVOLD --------
7fd59977 1683---- */
1684
1685 i__2 = ncf - 1;
1686 for (ncj = 1; ncj <= i__2; ++ncj) {
1687 bid = tabaux[ncj - 1];
1688 i__3 = *ndim;
1689 for (nd = 1; nd <= i__3; ++nd) {
1690 crvnew[nd + ncj * crvnew_dim1] += crvold[ncf + nd *
1691 crvold_dim1] * bid;
1692/* L400: */
1693 }
1694/* L300: */
1695 }
1696
1697 bid = tabaux[ncf - 1];
1698 i__2 = *ndim;
1699 for (nd = 1; nd <= i__2; ++nd) {
1700 crvnew[nd + ncf * crvnew_dim1] = crvold[ncf + nd * crvold_dim1] *
1701 bid;
1702/* L500: */
1703 }
1704
0d969553 1705/* --------- Calculation of (NCF+1) coeff. of [X1*t + X0]**(NCF) --------
7fd59977 1706---- */
1707
1708 tabaux[ncf] = tabaux[ncf - 1] * x1;
1709 for (ncj = ncf; ncj >= 2; --ncj) {
1710 tabaux[ncj - 1] = tabaux[ncj - 1] * x0 + tabaux[ncj - 2] * x1;
1711/* L600: */
1712 }
1713 tabaux[0] *= x0;
1714
1715/* L200: */
1716 }
1717
0d969553 1718/* -------------- Take into account the last coeff. of CRVOLD -----------
7fd59977 1719*/
1720
1721 i__1 = *ncoeff - 1;
1722 for (ncj = 1; ncj <= i__1; ++ncj) {
1723 bid = tabaux[ncj - 1];
1724 i__2 = *ndim;
1725 for (nd = 1; nd <= i__2; ++nd) {
1726 crvnew[nd + ncj * crvnew_dim1] += crvold[*ncoeff + nd *
1727 crvold_dim1] * bid;
1728/* L800: */
1729 }
1730/* L700: */
1731 }
1732 i__1 = *ndim;
1733 for (nd = 1; nd <= i__1; ++nd) {
1734 crvnew[nd + *ncoeff * crvnew_dim1] = crvold[*ncoeff + nd *
1735 crvold_dim1] * tabaux[*ncoeff - 1];
1736/* L900: */
1737 }
1738
1739/* ---------------------------- The end ---------------------------------
1740*/
1741
1742L9999:
1743 if (*iercod > 0) {
1744 AdvApp2Var_SysBase::maermsg_("MMARCIN", iercod, 7L);
1745 }
1746 if (ibb >= 2) {
1747 AdvApp2Var_SysBase::mgsomsg_("MMARCIN", 7L);
1748 }
1749 return 0;
1750} /* mmarcin_ */
1751
1752//=======================================================================
1753//function : mmatvec_
1754//purpose :
1755//=======================================================================
1756int mmatvec_(integer *nligne,
1757 integer *,//ncolon,
1758 integer *gposit,
1759 integer *,//gnstoc,
1760 doublereal *gmatri,
1761 doublereal *vecin,
1762 integer *deblig,
1763 doublereal *vecout,
1764 integer *iercod)
1765
1766{
1767 /* System generated locals */
1768 integer i__1, i__2;
1769
1770 /* Local variables */
1ef32e96
RL
1771 logical ldbg;
1772 integer jmin, jmax, i__, j, k;
1773 doublereal somme;
1774 integer aux;
7fd59977 1775
1776
1777/* ***********************************************************************
1778 */
1779
0d969553 1780/* FUNCTION : */
7fd59977 1781/* ---------- */
0d969553 1782/* Produce vector matrix in form of profile */
7fd59977 1783
1784
1785/* MOTS CLES : */
1786/* ----------- */
0d969553 1787/* RESERVE, MATRIX, PRODUCT, VECTOR, PROFILE */
7fd59977 1788
0d969553 1789/* INPUT ARGUMENTS : */
7fd59977 1790/* -------------------- */
0d969553
Y
1791/* NLIGNE : Line number of the matrix of constraints */
1792/* NCOLON : Number of column of the matrix of constraints */
1793/* GNSTOC: Number of coefficients in the profile of matrix GMATRI */
1794
1795/* GPOSIT: Table of positioning of terms of storage */
1796/* GPOSIT(1,I) contains the number of terms-1 on the line I
1797/* in the profile of the matrix. */
1798/* GPOSIT(2,I) contains the index of storage of diagonal term*/
1799/* of line I */
1800/* GPOSIT(3,I) contains the index of column of the first term of */
1801/* profile of line I */
1802/* GNSTOC: Number of coefficients in the profile of matrix */
7fd59977 1803/* GMATRI */
0d969553
Y
1804/* GMATRI : Matrix of constraints in form of profile */
1805/* VECIN : Input vector */
1806/* DEBLIG : Line indexusing which the vector matrix is calculated */
1807/*
1808/* OUTPUT ARGUMENTS */
7fd59977 1809/* --------------------- */
0d969553 1810/* VECOUT : VECTOR PRODUCT */
7fd59977 1811
0d969553 1812/* IERCOD : ERROR CODE */
7fd59977 1813
1814
0d969553 1815/* COMMONS USED : */
7fd59977 1816/* ------------------ */
1817
1818
0d969553 1819/* REFERENCES CALLED : */
7fd59977 1820/* --------------------- */
1821
1822
0d969553 1823/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 1824/* ----------------------------------- */
1825
7fd59977 1826/* ***********************************************************************
1827 */
1828/* DECLARATIONS */
1829/* ***********************************************************************
1830 */
1831
1832
1833
1834/* ***********************************************************************
1835 */
1836/* INITIALISATIONS */
1837/* ***********************************************************************
1838 */
1839
1840 /* Parameter adjustments */
1841 --vecout;
1842 gposit -= 4;
1843 --vecin;
1844 --gmatri;
1845
1846 /* Function Body */
1847 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1848 if (ldbg) {
1849 AdvApp2Var_SysBase::mgenmsg_("MMATVEC", 7L);
1850 }
1851 *iercod = 0;
1852
1853/* ***********************************************************************
1854 */
0d969553 1855/* Processing */
7fd59977 1856/* ***********************************************************************
1857 */
fadcea2c
RL
1858 AdvApp2Var_SysBase::mvriraz_(nligne,
1859 &vecout[1]);
7fd59977 1860 i__1 = *nligne;
1861 for (i__ = *deblig; i__ <= i__1; ++i__) {
1862 somme = 0.;
1863 jmin = gposit[i__ * 3 + 3];
1864 jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
1865 aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
1866 i__2 = jmax;
1867 for (j = jmin; j <= i__2; ++j) {
1868 k = j + aux;
1869 somme += gmatri[k] * vecin[j];
1870 }
1871 vecout[i__] = somme;
1872 }
1873
1874
1875
1876
1877
1878 goto L9999;
1879
1880/* ***********************************************************************
1881 */
0d969553 1882/* ERROR PROCESSING */
7fd59977 1883/* ***********************************************************************
1884 */
1885
1886
1887
1888
1889/* ***********************************************************************
1890 */
0d969553 1891/* RETURN CALLING PROGRAM */
7fd59977 1892/* ***********************************************************************
1893 */
1894
1895L9999:
1896
1897/* ___ DESALLOCATION, ... */
1898
1899 AdvApp2Var_SysBase::maermsg_("MMATVEC", iercod, 7L);
1900 if (ldbg) {
1901 AdvApp2Var_SysBase::mgsomsg_("MMATVEC", 7L);
1902 }
1903
1904 return 0 ;
1905} /* mmatvec_ */
1906
1907//=======================================================================
1908//function : mmbulld_
1909//purpose :
1910//=======================================================================
1911int AdvApp2Var_MathBase::mmbulld_(integer *nbcoln,
1912 integer *nblign,
1913 doublereal *dtabtr,
1914 integer *numcle)
1915
1916{
1917 /* System generated locals */
1918 integer dtabtr_dim1, dtabtr_offset, i__1, i__2;
1919
1920 /* Local variables */
1ef32e96
RL
1921 logical ldbg;
1922 doublereal daux;
1923 integer nite1, nite2, nchan, i1, i2;
7fd59977 1924
1925/* ***********************************************************************
1926 */
1927
0d969553 1928/* FUNCTION : */
7fd59977 1929/* ---------- */
0d969553
Y
1930/* Parsing of columns of a table of integers in increasing order */
1931/* KEYWORDS : */
7fd59977 1932/* ----------- */
0d969553
Y
1933/* POINT-ENTRY, PARSING */
1934/* INPUT ARGUMENTS : */
7fd59977 1935/* -------------------- */
0d969553
Y
1936/* - NBCOLN : Number of columns in the table */
1937/* - NBLIGN : Number of lines in the table */
1938/* - DTABTR : Table of integers to be parsed */
1939/* - NUMCLE : Position of the key on the column */
7fd59977 1940
0d969553 1941/* OUTPUT ARGUMENTS : */
7fd59977 1942/* --------------------- */
0d969553 1943/* - DTABTR : Parsed table */
7fd59977 1944
0d969553 1945/* COMMONS USED : */
7fd59977 1946/* ------------------ */
1947
1948
0d969553 1949/* REFERENCES CALLED : */
7fd59977 1950/* --------------------- */
1951
1952
0d969553 1953/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 1954/* ----------------------------------- */
0d969553
Y
1955/* Particularly performant if the table is almost parsed */
1956/* In the opposite case it is better to use MVSHELD */
7fd59977 1957/* ***********************************************************************
1958 */
1959
1960 /* Parameter adjustments */
1961 dtabtr_dim1 = *nblign;
1962 dtabtr_offset = dtabtr_dim1 + 1;
1963 dtabtr -= dtabtr_offset;
1964
1965 /* Function Body */
1966 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1967 if (ldbg) {
1968 AdvApp2Var_SysBase::mgenmsg_("MMBULLD", 7L);
1969 }
1970 nchan = 1;
1971 nite1 = *nbcoln;
1972 nite2 = 2;
1973
1974/* ***********************************************************************
1975 */
0d969553 1976/* PROCESSING */
7fd59977 1977/* ***********************************************************************
1978 */
1979
0d969553 1980/* ---->ALGORITHM in N^2 / 2 additional iteration */
7fd59977 1981
1982 while(nchan != 0) {
1983
0d969553 1984/* ----> Parsing from left to the right */
7fd59977 1985
1986 nchan = 0;
1987 i__1 = nite1;
1988 for (i1 = nite2; i1 <= i__1; ++i1) {
1989 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1 - 1)
1990 * dtabtr_dim1]) {
1991 i__2 = *nblign;
1992 for (i2 = 1; i2 <= i__2; ++i2) {
1993 daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
1994 dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
1995 dtabtr_dim1];
1996 dtabtr[i2 + i1 * dtabtr_dim1] = daux;
1997 }
1998 if (nchan == 0) {
1999 nchan = 1;
2000 }
2001 }
2002 }
2003 --nite1;
2004
0d969553 2005/* ----> Parsing from right to the left */
7fd59977 2006
2007 if (nchan != 0) {
2008 nchan = 0;
2009 i__1 = nite2;
2010 for (i1 = nite1; i1 >= i__1; --i1) {
2011 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1
2012 - 1) * dtabtr_dim1]) {
2013 i__2 = *nblign;
2014 for (i2 = 1; i2 <= i__2; ++i2) {
2015 daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
2016 dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
2017 dtabtr_dim1];
2018 dtabtr[i2 + i1 * dtabtr_dim1] = daux;
2019 }
2020 if (nchan == 0) {
2021 nchan = 1;
2022 }
2023 }
2024 }
2025 ++nite2;
2026 }
2027 }
2028
2029
2030 goto L9999;
2031
2032/* ***********************************************************************
2033 */
0d969553 2034/* ERROR PROCESSING */
7fd59977 2035/* ***********************************************************************
2036 */
2037
0d969553 2038/* ----> No errors at calling functions, only tests and loops. */
7fd59977 2039
2040/* ***********************************************************************
2041 */
0d969553 2042/* RETURN CALLING PROGRAM */
7fd59977 2043/* ***********************************************************************
2044 */
2045
2046L9999:
2047
2048 if (ldbg) {
2049 AdvApp2Var_SysBase::mgsomsg_("MMBULLD", 7L);
2050 }
2051
2052 return 0 ;
2053} /* mmbulld_ */
2054
2055
2056//=======================================================================
2057//function : AdvApp2Var_MathBase::mmcdriv_
2058//purpose :
2059//=======================================================================
2060int AdvApp2Var_MathBase::mmcdriv_(integer *ndimen,
2061 integer *ncoeff,
2062 doublereal *courbe,
2063 integer *ideriv,
2064 integer *ncofdv,
2065 doublereal *crvdrv)
2066
2067
2068{
2069 /* System generated locals */
2070 integer courbe_dim1, courbe_offset, crvdrv_dim1, crvdrv_offset, i__1,
2071 i__2;
2072
2073 /* Local variables */
1ef32e96
RL
2074 integer i__, j, k;
2075 doublereal mfactk, bid;
7fd59977 2076
2077
2078/* ***********************************************************************
2079 */
2080
0d969553 2081/* FUNCTION : */
7fd59977 2082/* ---------- */
0d969553
Y
2083/* Calculate matrix of a derivate curve of order IDERIV. */
2084/* with input parameters other than output parameters. */
7fd59977 2085
2086
0d969553 2087/* KEYWORDS : */
7fd59977 2088/* ----------- */
0d969553 2089/* COEFFICIENTS,CURVE,DERIVATE I-EME. */
7fd59977 2090
0d969553 2091/* INPUT ARGUMENTS : */
7fd59977 2092/* ------------------ */
0d969553
Y
2093/* NDIMEN : Space dimension (2 or 3 in general) */
2094/* NCOEFF : Degree +1 of the curve. */
2095/* COURBE : Table of coefficients of the curve. */
2096/* IDERIV : Required order of derivation : 1=1st derivate, etc... */
7fd59977 2097
0d969553 2098/* OUTPUT ARGUMENTS : */
7fd59977 2099/* ------------------- */
0d969553
Y
2100/* NCOFDV : Degree +1 of the derivative of order IDERIV of the curve. */
2101/* CRVDRV : Table of coefficients of the derivative of order IDERIV */
2102/* of the curve. */
7fd59977 2103
0d969553 2104/* COMMONS USED : */
7fd59977 2105/* ---------------- */
2106
0d969553 2107/* REFERENCES CALLED : */
7fd59977 2108/* ----------------------- */
2109
0d969553 2110/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 2111/* ----------------------------------- */
2112
0d969553
Y
2113/* ---> It is possible to take as output argument the curve */
2114/* and the number of coeff passed at input by making : */
7fd59977 2115/* CALL MMCDRIV(NDIMEN,NCOEFF,COURBE,IDERIV,NCOEFF,COURBE). */
0d969553
Y
2116/* After this call, NCOEFF does the number of coeff of the derived */
2117/* curve the coefficients which of are stored in CURVE. */
2118/* Attention to the coefficients of CURVE of rank superior to */
2119/* NCOEFF : they are not set to zero. */
7fd59977 2120
0d969553
Y
2121/* ---> Algorithm : */
2122/* The code below was written basing on the following algorithm:
7fd59977 2123*/
2124
0d969553
Y
2125/* Let P(t) = a1 + a2*t + ... an*t**n. Derivate of order k of P */
2126/* (containing n-k coefficients) is calculated as follows : */
7fd59977 2127
2128/* Pk(t) = a(k+1)*CNP(k,k)*k! */
2129/* + a(k+2)*CNP(k+1,k)*k! * t */
2130/* . */
2131/* . */
2132/* . */
2133/* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
7fd59977 2134/* ***********************************************************************
2135 */
2136
2137
0d969553 2138/* -------------- Case when the order of derivative is -------------------
7fd59977 2139*/
0d969553 2140/* ---------------- greater than the degree of the curve ---------------------
7fd59977 2141*/
2142
2143/* **********************************************************************
2144*/
2145
0d969553 2146/* FUNCTION : */
7fd59977 2147/* ---------- */
0d969553 2148/* Serves to provide the coefficients of binome (Pascal's triangle). */
7fd59977 2149
0d969553 2150/* KEYWORDS : */
7fd59977 2151/* ----------- */
0d969553 2152/* Binomial coeff from 0 to 60. read only . init par block data */
7fd59977 2153
0d969553 2154/* DEMSCRIPTION/NOTES/LIMITATIONS : */
7fd59977 2155/* ----------------------------------- */
0d969553
Y
2156/* Binomial coefficients form a triangular matrix. */
2157/* This matrix is completed in table CNP by its transposition. */
2158/* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
7fd59977 2159
0d969553
Y
2160/* Initialization is done by block-data MMLLL09.RES, */
2161/* created by program MQINICNP.FOR). */
7fd59977 2162/* **********************************************************************
2163*/
2164
2165
2166
2167/* ***********************************************************************
2168 */
2169
2170 /* Parameter adjustments */
2171 crvdrv_dim1 = *ndimen;
2172 crvdrv_offset = crvdrv_dim1 + 1;
2173 crvdrv -= crvdrv_offset;
2174 courbe_dim1 = *ndimen;
2175 courbe_offset = courbe_dim1 + 1;
2176 courbe -= courbe_offset;
2177
2178 /* Function Body */
2179 if (*ideriv >= *ncoeff) {
2180 i__1 = *ndimen;
2181 for (i__ = 1; i__ <= i__1; ++i__) {
2182 crvdrv[i__ + crvdrv_dim1] = 0.;
2183/* L10: */
2184 }
2185 *ncofdv = 1;
2186 goto L9999;
2187 }
2188/* **********************************************************************
2189*/
0d969553 2190/* General processing */
7fd59977 2191/* **********************************************************************
2192*/
0d969553 2193/* --------------------- Calculation of Factorial(IDERIV) ------------------
7fd59977 2194*/
2195
2196 k = *ideriv;
2197 mfactk = 1.;
2198 i__1 = k;
2199 for (i__ = 2; i__ <= i__1; ++i__) {
2200 mfactk *= i__;
2201/* L50: */
2202 }
2203
0d969553 2204/* ------------ Calculation of coeff of the derived of order IDERIV ----------
7fd59977 2205*/
0d969553
Y
2206/* ---> Attention : coefficient binomial C(n,m) is represented in */
2207/* MCCNP by CNP(N+1,M+1). */
7fd59977 2208
2209 i__1 = *ncoeff;
2210 for (j = k + 1; j <= i__1; ++j) {
2211 bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
2212 i__2 = *ndimen;
2213 for (i__ = 1; i__ <= i__2; ++i__) {
2214 crvdrv[i__ + (j - k) * crvdrv_dim1] = bid * courbe[i__ + j *
2215 courbe_dim1];
2216/* L200: */
2217 }
2218/* L100: */
2219 }
2220
2221 *ncofdv = *ncoeff - *ideriv;
2222
2223/* -------------------------------- The end -----------------------------
2224*/
2225
2226L9999:
2227 return 0;
2228} /* mmcdriv_ */
2229
2230//=======================================================================
2231//function : AdvApp2Var_MathBase::mmcglc1_
2232//purpose :
2233//=======================================================================
2234int AdvApp2Var_MathBase::mmcglc1_(integer *ndimax,
2235 integer *ndimen,
2236 integer *ncoeff,
2237 doublereal *courbe,
2238 doublereal *tdebut,
2239 doublereal *tfinal,
2240 doublereal *epsiln,
2241 doublereal *xlongc,
2242 doublereal *erreur,
2243 integer *iercod)
2244
2245
2246{
2247 /* System generated locals */
2248 integer courbe_dim1, courbe_offset, i__1;
2249 doublereal d__1;
2250
2251 /* Local variables */
1ef32e96
RL
2252 integer ndec;
2253 doublereal tdeb, tfin;
2254 integer iter;
2255 doublereal oldso;
2256 integer itmax;
2257 doublereal sottc;
2258 integer kk, ibb;
2259 doublereal dif, pas;
2260 doublereal som;
7fd59977 2261
2262
2263/* ***********************************************************************
2264 */
2265
0d969553 2266/* FUNCTION : */
7fd59977 2267/* ---------- */
0d969553
Y
2268/* Allows calculating the length of an arc of curve POLYNOMIAL */
2269/* on an interval [A,B]. */
7fd59977 2270
0d969553 2271/* KEYWORDS : */
7fd59977 2272/* ----------- */
0d969553 2273/* LENGTH,CURVE,GAUSS,PRIVATE. */
7fd59977 2274
0d969553 2275/* INPUT ARGUMENTS : */
7fd59977 2276/* ------------------ */
0d969553
Y
2277/* NDIMAX : Max. number of lines of tables */
2278/* (i.e. max. nb of polynoms). */
2279/* NDIMEN : Dimension of the space (nb of polynoms). */
2280/* NCOEFF : Nb of coefficients of the polynom. This is degree + 1.
2281*/
2282/* COURBE(NDIMAX,NCOEFF) : Coefficients of the curve. */
2283/* TDEBUT : Lower limit of the interval of integration for */
2284/* length calculation. */
2285/* TFINAL : Upper limit of the interval of integration for */
2286/* length calculation. */
2287/* EPSILN : REQIRED precision for length calculation. */
2288
2289/* OUTPUT ARGUMENTS : */
7fd59977 2290/* ------------------- */
0d969553
Y
2291/* XLONGC : Length of the arc of curve */
2292/* ERREUR : Precision OBTAINED for the length calculation. */
2293/* IERCOD : Error code, 0 OK, >0 Serious error. */
2294/* = 1 Too much iterations, the best calculated resultat */
2295/* (is almost ERROR) */
2296/* = 2 Pb MMLONCV (no result) */
2297/* = 3 NDIM or NCOEFF invalid (no result) */
2298
2299/* COMMONS USED : */
7fd59977 2300/* ---------------- */
2301
0d969553 2302/* REFERENCES CALLED : */
7fd59977 2303/* ----------------------- */
2304
0d969553 2305/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 2306/* ----------------------------------- */
0d969553
Y
2307/* The polynom is actually a set of polynoms with */
2308/* coefficients arranged in a table of 2 indices, */
2309/* each line relative to the polynom. */
2310/* The polynom is defined by these coefficients ordered */
2311/* by increasing power of the variable. */
2312/* All polynoms have the same number of coefficients (the */
2313/* same degree). */
2314
2315/* This program cancels and replaces LENGCV, MLONGC and MLENCV. */
2316
2317/* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
2318
7fd59977 2319/* > */
2320/* ***********************************************************************
2321 */
2322
0d969553 2323/* Name of the routine */
7fd59977 2324
2325
0d969553 2326/* ------------------------ General Initialization ---------------------
7fd59977 2327*/
2328
2329 /* Parameter adjustments */
2330 courbe_dim1 = *ndimax;
2331 courbe_offset = courbe_dim1 + 1;
2332 courbe -= courbe_offset;
2333
2334 /* Function Body */
2335 ibb = AdvApp2Var_SysBase::mnfndeb_();
2336 if (ibb >= 2) {
2337 AdvApp2Var_SysBase::mgenmsg_("MMCGLC1", 7L);
2338 }
2339
2340 *iercod = 0;
2341 *xlongc = 0.;
2342 *erreur = 0.;
2343
0d969553 2344/* ------ Test of equity of limits */
7fd59977 2345
2346 if (*tdebut == *tfinal) {
2347 *iercod = 0;
2348 goto L9999;
2349 }
2350
0d969553 2351/* ------ Test of the dimension and the number of coefficients */
7fd59977 2352
2353 if (*ndimen <= 0 || *ncoeff <= 0) {
2354 goto L9003;
2355 }
2356
0d969553
Y
2357/* ----- Nb of current cutting, nb of iteration, */
2358/* max nb of iterations */
7fd59977 2359
2360 ndec = 1;
2361 iter = 1;
2362
7fd59977 2363 itmax = 13;
2364
0d969553
Y
2365/* ------ Variation of the nb of intervals */
2366/* Multiplied by 2 at each iteration */
7fd59977 2367
2368L5000:
2369 pas = (*tfinal - *tdebut) / ndec;
2370 sottc = 0.;
2371
0d969553 2372/* ------ Loop on all current NDEC intervals */
7fd59977 2373
2374 i__1 = ndec;
2375 for (kk = 1; kk <= i__1; ++kk) {
2376
0d969553 2377/* ------ Limits of the current integration interval */
7fd59977 2378
2379 tdeb = *tdebut + (kk - 1) * pas;
2380 tfin = tdeb + pas;
2381 mmloncv_(ndimax, ndimen, ncoeff, &courbe[courbe_offset], &tdeb, &tfin,
2382 &som, iercod);
2383 if (*iercod > 0) {
2384 goto L9002;
2385 }
2386
2387 sottc += som;
2388
2389/* L100: */
2390 }
2391
2392
0d969553 2393/* ----------------- Test of the maximum number of iterations ------------
7fd59977 2394*/
2395
0d969553 2396/* Test if passes at least once ** */
7fd59977 2397
2398 if (iter == 1) {
2399 oldso = sottc;
2400 ndec <<= 1;
2401 ++iter;
2402 goto L5000;
2403 } else {
2404
0d969553 2405/* ------ Take into account DIF - Test of convergence */
7fd59977 2406
2407 ++iter;
41194117 2408 dif = (d__1 = sottc - oldso, advapp_abs(d__1));
7fd59977 2409
0d969553 2410/* ------ If DIF is OK, leave..., otherwise: */
7fd59977 2411
2412 if (dif > *epsiln) {
2413
0d969553 2414/* ------ If nb iteration exceeded, leave */
7fd59977 2415
2416 if (iter > itmax) {
2417 *iercod = 1;
2418 goto L9000;
2419 } else {
2420
0d969553 2421/* ------ Otherwise continue by cutting the initial interval.
7fd59977 2422 */
2423
2424 oldso = sottc;
2425 ndec <<= 1;
2426 goto L5000;
2427 }
2428 }
2429 }
2430
2431/* ------------------------------ THE END -------------------------------
2432*/
2433
2434L9000:
2435 *xlongc = sottc;
2436 *erreur = dif;
2437 goto L9999;
2438
0d969553 2439/* ---> PB in MMLONCV */
7fd59977 2440
2441L9002:
2442 *iercod = 2;
2443 goto L9999;
2444
0d969553 2445/* ---> NCOEFF or NDIM invalid. */
7fd59977 2446
2447L9003:
2448 *iercod = 3;
2449 goto L9999;
2450
2451L9999:
2452 if (*iercod > 0) {
2453 AdvApp2Var_SysBase::maermsg_("MMCGLC1", iercod, 7L);
2454 }
2455 if (ibb >= 2) {
2456 AdvApp2Var_SysBase::mgsomsg_("MMCGLC1", 7L);
2457 }
2458 return 0;
2459} /* mmcglc1_ */
2460
2461//=======================================================================
2462//function : mmchole_
2463//purpose :
2464//=======================================================================
2465int mmchole_(integer *,//mxcoef,
2466 integer *dimens,
2467 doublereal *amatri,
2468 integer *aposit,
2469 integer *posuiv,
2470 doublereal *chomat,
2471 integer *iercod)
2472
2473{
2474 /* System generated locals */
2475 integer i__1, i__2, i__3;
2476 doublereal d__1;
2477
2478 /* Builtin functions */
2479 //double sqrt();
2480
2481 /* Local variables */
1ef32e96
RL
2482 logical ldbg;
2483 integer kmin, i__, j, k;
2484 doublereal somme;
2485 integer ptini, ptcou;
7fd59977 2486
2487
2488/* ***********************************************************************
2489 */
2490
0d969553 2491/* FUNCTION : */
7fd59977 2492/* ---------- T */
0d969553
Y
2493/* Produce decomposition of choleski of matrix A in S.S */
2494/* Calculate inferior triangular matrix S. */
7fd59977 2495
0d969553 2496/* KEYWORDS : */
7fd59977 2497/* ----------- */
0d969553 2498/* RESOLUTION, MFACTORISATION, MATRIX_PROFILE, CHOLESKI */
7fd59977 2499
0d969553 2500/* INPUT ARGUMENTS : */
7fd59977 2501/* -------------------- */
0d969553
Y
2502/* MXCOEF : Max number of terms in the hessian profile */
2503/* DIMENS : Dimension of the problem */
2504/* AMATRI(MXCOEF) : Coefficients of the matrix profile */
2505/* APOSIT(1,*) : Distance diagonal-left extremity of the line
7fd59977 2506*/
0d969553
Y
2507/* APOSIT(2,*) : Position of diagonal terms in HESSIE */
2508/* POSUIV(MXCOEF) : first line inferior not out of profile */
7fd59977 2509
0d969553 2510/* OUTPUT ARGUMENTS : */
7fd59977 2511/* --------------------- */
0d969553
Y
2512/* CHOMAT(MXCOEF) : Inferior triangular matrix preserving the */
2513/* profile of AMATRI. */
2514/* IERCOD : error code */
7fd59977 2515/* = 0 : ok */
0d969553 2516/* = 1 : non-defined positive matrix */
7fd59977 2517
0d969553 2518/* COMMONS USED : */
7fd59977 2519/* ------------------ */
2520
2521/* .Neant. */
2522
0d969553 2523/* REFERENCES CALLED : */
7fd59977 2524/* ---------------------- */
2525
0d969553 2526/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 2527/* ----------------------------------- */
0d969553 2528/* DEBUG LEVEL = 4 */
7fd59977 2529/* ***********************************************************************
2530 */
2531/* DECLARATIONS */
2532/* ***********************************************************************
2533 */
2534
2535
2536
2537/* ***********************************************************************
2538 */
2539/* INITIALISATIONS */
2540/* ***********************************************************************
2541 */
2542
2543 /* Parameter adjustments */
2544 --chomat;
2545 --posuiv;
2546 --amatri;
2547 aposit -= 3;
2548
2549 /* Function Body */
2550 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
2551 if (ldbg) {
2552 AdvApp2Var_SysBase::mgenmsg_("MMCHOLE", 7L);
2553 }
2554 *iercod = 0;
2555
2556/* ***********************************************************************
2557 */
0d969553 2558/* PROCESSING */
7fd59977 2559/* ***********************************************************************
2560 */
2561
2562 i__1 = *dimens;
2563 for (j = 1; j <= i__1; ++j) {
2564
2565 ptini = aposit[(j << 1) + 2];
2566
2567 somme = 0.;
2568 i__2 = ptini - 1;
2569 for (k = ptini - aposit[(j << 1) + 1]; k <= i__2; ++k) {
2570/* Computing 2nd power */
2571 d__1 = chomat[k];
2572 somme += d__1 * d__1;
2573 }
2574
2575 if (amatri[ptini] - somme < 1e-32) {
2576 goto L9101;
2577 }
2578 chomat[ptini] = sqrt(amatri[ptini] - somme);
2579
2580 ptcou = ptini;
2581
2582 while(posuiv[ptcou] > 0) {
2583
2584 i__ = posuiv[ptcou];
2585 ptcou = aposit[(i__ << 1) + 2] - (i__ - j);
2586
0d969553 2587/* Calculate the sum of S .S for k =1 a j-1 */
7fd59977 2588/* ik jk */
2589 somme = 0.;
2590/* Computing MAX */
2591 i__2 = i__ - aposit[(i__ << 1) + 1], i__3 = j - aposit[(j << 1) +
2592 1];
41194117 2593 kmin = advapp_max(i__2,i__3);
7fd59977 2594 i__2 = j - 1;
2595 for (k = kmin; k <= i__2; ++k) {
2596 somme += chomat[aposit[(i__ << 1) + 2] - (i__ - k)] * chomat[
2597 aposit[(j << 1) + 2] - (j - k)];
2598 }
2599
2600 chomat[ptcou] = (amatri[ptcou] - somme) / chomat[ptini];
2601 }
2602 }
2603
2604 goto L9999;
2605
2606/* ***********************************************************************
2607 */
0d969553 2608/* ERROR PROCESSING */
7fd59977 2609/* ***********************************************************************
2610 */
2611
2612L9101:
2613 *iercod = 1;
2614 goto L9999;
2615
2616/* ***********************************************************************
2617 */
0d969553 2618/* RETURN CALLING PROGRAM */
7fd59977 2619/* ***********************************************************************
2620 */
2621
2622L9999:
2623
2624 AdvApp2Var_SysBase::maermsg_("MMCHOLE", iercod, 7L);
2625 if (ldbg) {
2626 AdvApp2Var_SysBase::mgsomsg_("MMCHOLE", 7L);
2627 }
2628
2629 return 0 ;
2630} /* mmchole_ */
2631
2632//=======================================================================
2633//function : AdvApp2Var_MathBase::mmcvctx_
2634//purpose :
2635//=======================================================================
2636int AdvApp2Var_MathBase::mmcvctx_(integer *ndimen,
2637 integer *ncofmx,
2638 integer *nderiv,
2639 doublereal *ctrtes,
2640 doublereal *crvres,
2641 doublereal *tabaux,
2642 doublereal *xmatri,
2643 integer *iercod)
2644
2645{
2646 /* System generated locals */
2647 integer ctrtes_dim1, ctrtes_offset, crvres_dim1, crvres_offset,
2648 xmatri_dim1, xmatri_offset, tabaux_dim1, tabaux_offset, i__1,
2649 i__2;
2650
2651 /* Local variables */
1ef32e96
RL
2652 integer moup1, nordr;
2653 integer nd;
2654 integer ibb, ncf, ndv;
2655 doublereal eps1;
7fd59977 2656
2657
2658/* ***********************************************************************
2659 */
2660
0d969553 2661/* FUNCTION : */
7fd59977 2662/* ---------- */
0d969553
Y
2663/* Calculate a polynomial curve checking the */
2664/* passage constraints (interpolation) */
2665/* from first derivatives, etc... to extremities. */
2666/* Parameters at the extremities are supposed to be -1 and 1. */
7fd59977 2667
0d969553 2668/* KEYWORDS : */
7fd59977 2669/* ----------- */
0d969553 2670/* ALL, AB_SPECIFI::CONSTRAINTS&,INTERPOLATION,&CURVE */
7fd59977 2671
0d969553 2672/* INPUT ARGUMENTS : */
7fd59977 2673/* ------------------ */
0d969553
Y
2674/* NDIMEN : Space Dimension. */
2675/* NCOFMX : Nb of coeff. of curve CRVRES on each */
7fd59977 2676/* dimension. */
0d969553 2677/* NDERIV : Order of constraint with derivatives : */
7fd59977 2678/* 0 --> interpolation simple. */
0d969553
Y
2679/* 1 --> interpolation+constraints with 1st. */
2680/* 2 --> cas (0)+ (1) + " " 2nd derivatives. */
7fd59977 2681/* etc... */
0d969553
Y
2682/* CTRTES : Table of constraints. */
2683/* CTRTES(*,1,*) = contraints at -1. */
2684/* CTRTES(*,2,*) = contraints at 1. */
7fd59977 2685
0d969553 2686/* OUTPUT ARGUMENTS : */
7fd59977 2687/* ------------------- */
0d969553
Y
2688/* CRVRES : Resulting curve defined on (-1,1). */
2689/* TABAUX : Auxilliary matrix. */
2690/* XMATRI : Auxilliary matrix. */
7fd59977 2691
2692/* COMMONS UTILISES : */
2693/* ---------------- */
2694
2695/* .Neant. */
2696
0d969553 2697/* REFERENCES CALLED : */
7fd59977 2698/* ---------------------- */
2699/* Type Name */
2700/* MAERMSG R*8 DFLOAT MGENMSG */
2701/* MGSOMSG MMEPS1 MMRSLW */
2702/* I*4 MNFNDEB */
2703
0d969553 2704/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 2705/* ----------------------------------- */
0d969553
Y
2706/* The polynom (or the curve) is calculated by solving a */
2707/* system of linear equations. If the imposed degree is great */
2708/* it is preferable to call a routine based on */
2709/* Lagrange or Hermite interpolation depending on the case. */
2710/* (for a high degree the matrix of the system can be badly */
2711/* conditionned). */
2712/* This routine returns a curve defined in (-1,1). */
2713/* In general case, it is necessary to use MCVCTG. */
7fd59977 2714/* > */
2715/* ***********************************************************************
2716 */
2717
0d969553 2718/* Name of the routine */
7fd59977 2719
2720
2721 /* Parameter adjustments */
2722 crvres_dim1 = *ncofmx;
2723 crvres_offset = crvres_dim1 + 1;
2724 crvres -= crvres_offset;
2725 xmatri_dim1 = *nderiv + 1;
2726 xmatri_offset = xmatri_dim1 + 1;
2727 xmatri -= xmatri_offset;
2728 tabaux_dim1 = *nderiv + 1 + *ndimen;
2729 tabaux_offset = tabaux_dim1 + 1;
2730 tabaux -= tabaux_offset;
2731 ctrtes_dim1 = *ndimen;
2732 ctrtes_offset = ctrtes_dim1 * 3 + 1;
2733 ctrtes -= ctrtes_offset;
2734
2735 /* Function Body */
2736 ibb = AdvApp2Var_SysBase::mnfndeb_();
2737 if (ibb >= 3) {
2738 AdvApp2Var_SysBase::mgenmsg_("MMCVCTX", 7L);
2739 }
0d969553 2740/* Precision. */
7fd59977 2741 AdvApp2Var_MathBase::mmeps1_(&eps1);
2742
0d969553 2743/* ****************** CALCULATION OF EVEN COEFFICIENTS *********************
7fd59977 2744*/
0d969553 2745/* ------------------------- Initialization -----------------------------
7fd59977 2746*/
2747
2748 nordr = *nderiv + 1;
2749 i__1 = nordr;
2750 for (ncf = 1; ncf <= i__1; ++ncf) {
2751 tabaux[ncf + tabaux_dim1] = 1.;
2752/* L100: */
2753 }
2754
0d969553 2755/* ---------------- Calculation of terms corresponding to derivatives -------
7fd59977 2756*/
2757
2758 i__1 = nordr;
2759 for (ndv = 2; ndv <= i__1; ++ndv) {
2760 i__2 = nordr;
2761 for (ncf = 1; ncf <= i__2; ++ncf) {
2762 tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
2763 tabaux_dim1] * (doublereal) ((ncf << 1) - ndv);
2764/* L300: */
2765 }
2766/* L200: */
2767 }
2768
0d969553 2769/* ------------------ Writing the second member -----------------------
7fd59977 2770*/
2771
2772 moup1 = 1;
2773 i__1 = nordr;
2774 for (ndv = 1; ndv <= i__1; ++ndv) {
2775 i__2 = *ndimen;
2776 for (nd = 1; nd <= i__2; ++nd) {
2777 tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
2778 + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2779 * ctrtes_dim1]) / 2.;
2780/* L500: */
2781 }
2782 moup1 = -moup1;
2783/* L400: */
2784 }
2785
0d969553 2786/* -------------------- Resolution of the system ---------------------------
7fd59977 2787*/
2788
2789 mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2790 xmatri_offset], iercod);
2791 if (*iercod > 0) {
2792 goto L9999;
2793 }
2794 i__1 = *ndimen;
2795 for (nd = 1; nd <= i__1; ++nd) {
2796 i__2 = nordr;
2797 for (ncf = 1; ncf <= i__2; ++ncf) {
2798 crvres[(ncf << 1) - 1 + nd * crvres_dim1] = xmatri[ncf + nd *
2799 xmatri_dim1];
2800/* L700: */
2801 }
2802/* L600: */
2803 }
2804
0d969553 2805/* ***************** CALCULATION OF UNEVEN COEFFICIENTS ********************
7fd59977 2806*/
0d969553 2807/* ------------------------- Initialization -----------------------------
7fd59977 2808*/
2809
2810
2811 i__1 = nordr;
2812 for (ncf = 1; ncf <= i__1; ++ncf) {
2813 tabaux[ncf + tabaux_dim1] = 1.;
2814/* L1100: */
2815 }
2816
0d969553 2817/* ---------------- Calculation of terms corresponding to derivatives -------
7fd59977 2818*/
2819
2820 i__1 = nordr;
2821 for (ndv = 2; ndv <= i__1; ++ndv) {
2822 i__2 = nordr;
2823 for (ncf = 1; ncf <= i__2; ++ncf) {
2824 tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
2825 tabaux_dim1] * (doublereal) ((ncf << 1) - ndv + 1);
2826/* L1300: */
2827 }
2828/* L1200: */
2829 }
2830
0d969553 2831/* ------------------ Writing of the second member -----------------------
7fd59977 2832*/
2833
2834 moup1 = -1;
2835 i__1 = nordr;
2836 for (ndv = 1; ndv <= i__1; ++ndv) {
2837 i__2 = *ndimen;
2838 for (nd = 1; nd <= i__2; ++nd) {
2839 tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
2840 + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2841 * ctrtes_dim1]) / 2.;
2842/* L1500: */
2843 }
2844 moup1 = -moup1;
2845/* L1400: */
2846 }
2847
0d969553 2848/* -------------------- Solution of the system ---------------------------
7fd59977 2849*/
2850
2851 mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2852 xmatri_offset], iercod);
2853 if (*iercod > 0) {
2854 goto L9999;
2855 }
2856 i__1 = *ndimen;
2857 for (nd = 1; nd <= i__1; ++nd) {
2858 i__2 = nordr;
2859 for (ncf = 1; ncf <= i__2; ++ncf) {
2860 crvres[(ncf << 1) + nd * crvres_dim1] = xmatri[ncf + nd *
2861 xmatri_dim1];
2862/* L1700: */
2863 }
2864/* L1600: */
2865 }
2866
2867/* --------------------------- The end ----------------------------------
2868*/
2869
2870L9999:
2871 if (*iercod != 0) {
2872 AdvApp2Var_SysBase::maermsg_("MMCVCTX", iercod, 7L);
2873 }
2874 if (ibb >= 3) {
2875 AdvApp2Var_SysBase::mgsomsg_("MMCVCTX", 7L);
2876 }
2877
2878 return 0 ;
2879} /* mmcvctx_ */
2880
2881//=======================================================================
2882//function : AdvApp2Var_MathBase::mmcvinv_
2883//purpose :
2884//=======================================================================
2885 int AdvApp2Var_MathBase::mmcvinv_(integer *ndimax,
2886 integer *ncoef,
2887 integer *ndim,
2888 doublereal *curveo,
2889 doublereal *curve)
2890
2891{
2892 /* Initialized data */
2893
2894 static char nomprg[8+1] = "MMCVINV ";
2895
2896 /* System generated locals */
2897 integer curve_dim1, curve_offset, curveo_dim1, curveo_offset, i__1, i__2;
2898
2899 /* Local variables */
1ef32e96 2900 integer i__, nd, ibb;
7fd59977 2901
2902
2903/* ***********************************************************************
2904 */
2905
0d969553 2906/* FUNCTION : */
7fd59977 2907/* ---------- */
0d969553 2908/* Inversion of arguments of the final curve. */
7fd59977 2909
0d969553 2910/* KEYWORDS : */
7fd59977 2911/* ----------- */
0d969553 2912/* SMOOTHING,CURVE */
7fd59977 2913
2914
0d969553 2915/* INPUT ARGUMENTS : */
7fd59977 2916/* ------------------ */
2917
0d969553
Y
2918/* NDIM: Space Dimension. */
2919/* NCOEF: Degree of the polynom. */
2920/* CURVEO: The curve before inversion. */
7fd59977 2921
0d969553 2922/* OUTPUT ARGUMENTS : */
7fd59977 2923/* ------------------- */
0d969553 2924/* CURVE: The curve after inversion. */
7fd59977 2925
0d969553 2926/* COMMONS USED : */
7fd59977 2927/* ---------------- */
7fd59977 2928/* REFERENCES APPELEES : */
2929/* ----------------------- */
0d969553 2930/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 2931/* ----------------------------------- */
7fd59977 2932/* ***********************************************************************
2933 */
2934
0d969553 2935/* The name of the routine */
7fd59977 2936 /* Parameter adjustments */
2937 curve_dim1 = *ndimax;
2938 curve_offset = curve_dim1 + 1;
2939 curve -= curve_offset;
2940 curveo_dim1 = *ncoef;
2941 curveo_offset = curveo_dim1 + 1;
2942 curveo -= curveo_offset;
2943
2944 /* Function Body */
2945
2946 ibb = AdvApp2Var_SysBase::mnfndeb_();
2947 if (ibb >= 2) {
2948 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
2949 }
2950
2951 i__1 = *ncoef;
2952 for (i__ = 1; i__ <= i__1; ++i__) {
2953 i__2 = *ndim;
2954 for (nd = 1; nd <= i__2; ++nd) {
2955 curve[nd + i__ * curve_dim1] = curveo[i__ + nd * curveo_dim1];
2956/* L300: */
2957 }
2958 }
2959
2960/* L9999: */
2961 return 0;
2962} /* mmcvinv_ */
2963
2964//=======================================================================
2965//function : mmcvstd_
2966//purpose :
2967//=======================================================================
2968int mmcvstd_(integer *ncofmx,
2969 integer *ndimax,
2970 integer *ncoeff,
2971 integer *ndimen,
2972 doublereal *crvcan,
2973 doublereal *courbe)
2974
2975{
2976 /* System generated locals */
2977 integer courbe_dim1, crvcan_dim1, crvcan_offset, i__1, i__2, i__3;
2978
2979 /* Local variables */
1ef32e96
RL
2980 integer ndeg, i__, j, j1, nd, ibb;
2981 doublereal bid;
7fd59977 2982
2983
2984/* ***********************************************************************
2985 */
2986
0d969553 2987/* FUNCTION : */
7fd59977 2988/* ---------- */
0d969553 2989/* Transform curve defined between [-1,1] into [0,1]. */
7fd59977 2990
0d969553 2991/* KEYWORDS : */
7fd59977 2992/* ----------- */
0d969553 2993/* LIMITATION,RESTRICTION,CURVE */
7fd59977 2994
0d969553 2995/* INPUT ARGUMENTS : */
7fd59977 2996/* ------------------ */
0d969553
Y
2997/* NDIMAX : Dimension of the space. */
2998/* NDIMEN : Dimension of the curve. */
2999/* NCOEFF : Degree of the curve. */
3000/* CRVCAN(NCOFMX,NDIMEN): The curve is defined at the interval [-1,1]. */
7fd59977 3001
0d969553 3002/* OUTPUT ARGUMENTS : */
7fd59977 3003/* ------------------- */
0d969553 3004/* CURVE(NDIMAX,NCOEFF): Curve defined at the interval [0,1]. */
7fd59977 3005
0d969553 3006/* COMMONS USED : */
7fd59977 3007/* ---------------- */
3008
0d969553 3009/* REFERENCES CALLED : */
7fd59977 3010/* ----------------------- */
3011
0d969553 3012/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 3013/* ----------------------------------- */
7fd59977 3014/* > */
3015/* ***********************************************************************
3016 */
3017
0d969553 3018/* Name of the program. */
7fd59977 3019
3020
3021/* **********************************************************************
3022*/
3023
0d969553 3024/* FUNCTION : */
7fd59977 3025/* ---------- */
0d969553 3026/* Provides binomial coefficients (Pascal triangle). */
7fd59977 3027
0d969553 3028/* KEYWORDS : */
7fd59977 3029/* ----------- */
0d969553 3030/* Binomial coefficient from 0 to 60. read only . init by block data */
7fd59977 3031
0d969553 3032/* DEMSCRIPTION/NOTES/LIMITATIONS : */
7fd59977 3033/* ----------------------------------- */
0d969553
Y
3034/* Binomial coefficients form a triangular matrix. */
3035/* This matrix is completed in table CNP by its transposition. */
3036/* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
7fd59977 3037
0d969553
Y
3038/* Initialization is done with block-data MMLLL09.RES, */
3039/* created by the program MQINICNP.FOR. */
7fd59977 3040/* > */
3041/* **********************************************************************
3042*/
3043
3044
3045
3046/* ***********************************************************************
3047 */
3048
3049 /* Parameter adjustments */
3050 courbe_dim1 = *ndimax;
3051 --courbe;
3052 crvcan_dim1 = *ncofmx;
3053 crvcan_offset = crvcan_dim1;
3054 crvcan -= crvcan_offset;
3055
3056 /* Function Body */
3057 ibb = AdvApp2Var_SysBase::mnfndeb_();
3058 if (ibb >= 3) {
3059 AdvApp2Var_SysBase::mgenmsg_("MMCVSTD", 7L);
3060 }
3061 ndeg = *ncoeff - 1;
3062
0d969553 3063/* ------------------ Construction of the resulting curve ----------------
7fd59977 3064*/
3065
3066 i__1 = *ndimen;
3067 for (nd = 1; nd <= i__1; ++nd) {
3068 i__2 = ndeg;
3069 for (j = 0; j <= i__2; ++j) {
3070 bid = 0.;
3071 i__3 = ndeg;
3072 for (i__ = j; i__ <= i__3; i__ += 2) {
3073 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
3074 * 61];
3075/* L410: */
3076 }
3077 courbe[nd + j * courbe_dim1] = bid;
3078
3079 bid = 0.;
3080 j1 = j + 1;
3081 i__3 = ndeg;
3082 for (i__ = j1; i__ <= i__3; i__ += 2) {
3083 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
3084 * 61];
3085/* L420: */
3086 }
3087 courbe[nd + j * courbe_dim1] -= bid;
3088/* L400: */
3089 }
3090/* L300: */
3091 }
3092
0d969553 3093/* ------------------- Renormalization of the CURVE -------------------------
7fd59977 3094 */
3095
3096 bid = 1.;
3097 i__1 = ndeg;
3098 for (i__ = 0; i__ <= i__1; ++i__) {
3099 i__2 = *ndimen;
3100 for (nd = 1; nd <= i__2; ++nd) {
3101 courbe[nd + i__ * courbe_dim1] *= bid;
3102/* L510: */
3103 }
3104 bid *= 2.;
3105/* L500: */
3106 }
3107
3108/* ----------------------------- The end --------------------------------
3109*/
3110
3111 if (ibb >= 3) {
3112 AdvApp2Var_SysBase::mgsomsg_("MMCVSTD", 7L);
3113 }
3114 return 0;
3115} /* mmcvstd_ */
3116
3117//=======================================================================
3118//function : AdvApp2Var_MathBase::mmdrc11_
3119//purpose :
3120//=======================================================================
3121int AdvApp2Var_MathBase::mmdrc11_(integer *iordre,
3122 integer *ndimen,
3123 integer *ncoeff,
3124 doublereal *courbe,
3125 doublereal *points,
3126 doublereal *mfactab)
3127
3128{
3129 /* System generated locals */
3130 integer courbe_dim1, courbe_offset, points_dim2, points_offset, i__1,
3131 i__2;
3132
3133 /* Local variables */
3134
1ef32e96 3135 integer ndeg, i__, j, ndgcb, nd, ibb;
7fd59977 3136
3137
3138/* **********************************************************************
3139*/
3140
0d969553 3141/* FUNCTION : */
7fd59977 3142/* ---------- */
0d969553
Y
3143/* Calculation of successive derivatives of equation CURVE with */
3144/* parameters -1, 1 from order 0 to order IORDRE */
3145/* included. The calculation is produced without knowing the coefficients of */
3146/* derivatives of the curve. */
7fd59977 3147
0d969553 3148/* KEYWORDS : */
7fd59977 3149/* ----------- */
0d969553 3150/* POSITIONING,EXTREMITIES,CURVE,DERIVATIVE. */
7fd59977 3151
0d969553 3152/* INPUT ARGUMENTS : */
7fd59977 3153/* ------------------ */
0d969553
Y
3154/* IORDRE : Maximum order of calculation of derivatives. */
3155/* NDIMEN : Dimension of the space. */
3156/* NCOEFF : Number of coefficients of the curve (degree+1). */
3157/* COURBE : Table of coefficients of the curve. */
7fd59977 3158
0d969553 3159/* OUTPUT ARGUMENTS : */
7fd59977 3160/* ------------------- */
0d969553
Y
3161/* POINTS : Table of values of consecutive derivatives */
3162/* of parameters -1.D0 and 1.D0. */
3163/* MFACTAB : Auxiliary table for calculation of factorial(I).
7fd59977 3164*/
3165
0d969553 3166/* COMMONS USED : */
7fd59977 3167/* ---------------- */
0d969553 3168/* None. */
7fd59977 3169
0d969553 3170/* REFERENCES CALLED : */
7fd59977 3171/* ----------------------- */
3172
0d969553 3173/* DESCRIPTION/NOTES/LIMITATIONS : */
7fd59977 3174/* ----------------------------------- */
3175
0d969553
Y
3176/* ---> ATTENTION, the coefficients of the curve are */
3177/* in a reverse order. */
7fd59977 3178
0d969553
Y
3179/* ---> The algorithm of calculation of derivatives is based on */
3180/* generalization of Horner scheme : */
7fd59977 3181/* k 2 */
0d969553 3182/* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
7fd59977 3183
3184
0d969553 3185/* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
7fd59977 3186
3187/* aj = a(j-1).x + u(k-j) */
3188/* bj = b(j-1).x + a(j-1) */
3189/* cj = c(j-1).x + b(j-1) */
3190
0d969553 3191/* So : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
7fd59977 3192
0d969553 3193/* The algorithm is generalized easily for calculation of */
7fd59977 3194
3195/* (n) */
3196/* C (x) . */
3197/* --------- */
3198/* n! */
3199
3200/* Reference : D. KNUTH, "The Art of Computer Programming" */
3201/* --------- Vol. 2/Seminumerical Algorithms */
3202/* Addison-Wesley Pub. Co. (1969) */
3203/* pages 423-425. */
7fd59977 3204/* > */
3205/* **********************************************************************
3206*/
3207
0d969553 3208/* Name of the routine */
7fd59977 3209
3210 /* Parameter adjustments */
3211 points_dim2 = *iordre + 1;
3212 points_offset = (points_dim2 << 1) + 1;
3213 points -= points_offset;
3214 courbe_dim1 = *ncoeff;
3215 courbe_offset = courbe_dim1;
3216 courbe -= courbe_offset;
3217
3218 /* Function Body */
3219 ibb = AdvApp2Var_SysBase::mnfndeb_();
3220 if (ibb >= 2) {
3221 AdvApp2Var_SysBase::mgenmsg_("MMDRC11", 7L);
3222 }
3223
3224 if (*iordre < 0 || *ncoeff < 1) {
3225 goto L9999;
3226 }
3227
0d969553 3228/* ------------------- Initialization of table POINTS -----------------
7fd59977 3229*/
3230
3231 ndgcb = *ncoeff - 1;
3232 i__1 = *ndimen;
3233 for (nd = 1; nd <= i__1; ++nd) {
3234 points[(nd * points_dim2 << 1) + 1] = courbe[ndgcb + nd * courbe_dim1]
3235 ;
3236 points[(nd * points_dim2 << 1) + 2] = courbe[ndgcb + nd * courbe_dim1]
3237 ;
3238/* L100: */
3239 }
3240
3241 i__1 = *ndimen;
3242 for (nd = 1; nd <= i__1; ++nd) {
3243 i__2 = *iordre;
3244 for (j = 1; j <= i__2; ++j) {
3245 points[((j + nd * points_dim2) << 1) + 1] = 0.;
3246 points[((j + nd * points_dim2) << 1) + 2] = 0.;
3247/* L400: */
3248 }
3249/* L300: */
3250 }
3251
0d969553 3252/* Calculation with parameter -1 and 1 */
7fd59977 3253
3254 i__1 = *ndimen;
3255 for (nd = 1; nd <= i__1; ++nd) {
3256 i__2 = ndgcb;
3257 for (ndeg = 1; ndeg <= i__2; ++ndeg) {
3258 for (i__ = *iordre; i__ >= 1; --i__) {
3259 points[((i__ + nd * points_dim2) << 1) + 1] = -points[((i__ + nd
3260 * points_dim2) << 1) + 1] + points[((i__ - 1 + nd *
3261 points_dim2) << 1) + 1];
3262 points[((i__ + nd * points_dim2) << 1) + 2] += points[((i__ - 1
3263 + nd * points_dim2) << 1) + 2];
3264/* L800: */
3265 }
3266 points[(nd * points_dim2 << 1) + 1] = -points[(nd * points_dim2 <<
3267 1) + 1] + courbe[ndgcb - ndeg + nd * courbe_dim1];
3268 points[(nd * points_dim2 << 1) + 2] += courbe[ndgcb - ndeg + nd *
3269 courbe_dim1];
3270/* L700: */
3271 }
3272/* L600: */
3273 }
3274