1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 // This file is part of Open CASCADE Technology software library.
5 // This library is free software; you can redistribute it and / or modify it
6 // under the terms of the GNU Lesser General Public version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
14 // AdvApp2Var_MathBase.cxx
16 #include <AdvApp2Var_SysBase.hxx>
17 #include <AdvApp2Var_Data_f2c.hxx>
18 #include <AdvApp2Var_MathBase.hxx>
19 #include <AdvApp2Var_Data.hxx>
23 int mmchole_(integer *mxcoef,
35 int mmrslss_(integer *mxcoef,
45 int mfac_(doublereal *f,
49 int mmaper0_(integer *ncofmx,
57 int mmaper2_(integer *ncofmx,
66 int mmaper4_(integer *ncofmx,
75 int mmaper6_(integer *ncofmx,
84 int mmarc41_(integer *ndimax,
94 int mmatvec_(integer *nligne,
105 int mmcvstd_(integer *ncofmx,
113 int mmdrvcb_(integer *ideriv,
122 int mmexthi_(integer *ndegre,
126 int mmextrl_(integer *ndegre,
132 int mmherm0_(doublereal *debfin,
136 int mmherm1_(doublereal *debfin,
142 int mmloncv_(integer *ndimax,
151 int mmpojac_(doublereal *tparam,
159 int mmrslw_(integer *normax,
167 int mmtmave_(integer *nligne,
176 int mmtrpj0_(integer *ncofmx,
185 int mmtrpj2_(integer *ncofmx,
195 int mmtrpj4_(integer *ncofmx,
204 int mmtrpj6_(integer *ncofmx,
213 integer pow__ii(integer *x,
217 int mvcvin2_(integer *ncoeff,
223 int mvcvinv_(integer *ncoeff,
229 int mvgaus0_(integer *kindic,
235 int mvpscr2_(integer *ncoeff,
241 int mvpscr3_(integer *ncoeff,
247 doublereal eps1, eps2, eps3, eps4;
248 integer niterm, niterr;
252 doublereal tdebut, tfinal, verifi, cmherm[576];
255 //=======================================================================
256 //function : AdvApp2Var_MathBase::mdsptpt_
258 //=======================================================================
259 int AdvApp2Var_MathBase::mdsptpt_(integer *ndimen,
266 /* System generated locals */
270 /* Local variables */
272 doublereal* differ = 0;
276 /* **********************************************************************
281 /* CALCULATE DISTANCE BETWEEN TWO POINTS */
285 /* DISTANCE,POINT. */
287 /* INPUT ARGUMENTS : */
288 /* ------------------ */
289 /* NDIMEN: Space Dimension. */
290 /* POINT1: Table of coordinates of the 1st point. */
291 /* POINT2: Table of coordinates of the 2nd point. */
293 /* OUTPUT ARGUMENTS : */
294 /* ------------------- */
295 /* DISTAN: Distance between 2 points. */
298 /* ---------------- */
300 /* REFERENCES CALLED : */
301 /* ----------------------- */
303 /* DESCRIPTION/NOTES/LIMITATIONS : */
304 /* ----------------------------------- */
306 /* **********************************************************************
310 /* ***********************************************************************
313 /* ***********************************************************************
316 /* Parameter adjustment */
324 /* ***********************************************************************
327 /* ***********************************************************************
330 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
332 anAdvApp2Var_SysBase.mcrrqst_(&c__8, ndimen, differ, &iofset, &ier);
335 /* --- If allocation is refused, the trivial method is applied. */
341 for (i__ = 1; i__ <= i__1; ++i__) {
342 /* Computing 2nd power */
343 d__1 = point1[i__] - point2[i__];
344 *distan += d__1 * d__1;
346 *distan = sqrt(*distan);
348 /* --- Otherwise MZSNORM is used to minimize the risks of overflow
353 for (i__ = 1; i__ <= i__1; ++i__) {
355 differ[j] = point2[i__] - point1[i__];
358 *distan = AdvApp2Var_MathBase::mzsnorm_(ndimen, &differ[iofset]);
362 /* ***********************************************************************
364 /* RETURN CALLING PROGRAM */
365 /* ***********************************************************************
368 /* --- Dynamic Desallocation */
371 anAdvApp2Var_SysBase.mcrdelt_(&c__8, ndimen, differ, &iofset, &ier);
377 //=======================================================================
380 //=======================================================================
381 int mfac_(doublereal *f,
385 /* System generated locals */
388 /* Local variables */
391 /* FORTRAN CONFORME AU TEXT */
392 /* CALCUL DE MFACTORIEL N */
393 /* Parameter adjustments */
399 for (i__ = 2; i__ <= i__1; ++i__) {
401 f[i__] = i__ * f[i__ - 1];
406 //=======================================================================
407 //function : AdvApp2Var_MathBase::mmapcmp_
409 //=======================================================================
410 int AdvApp2Var_MathBase::mmapcmp_(integer *ndim,
417 /* System generated locals */
418 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
421 /* Local variables */
422 integer ipair, nd, ndegre, impair, ibb, idg;
423 //extern int mgsomsg_();//mgenmsg_(),
425 /* **********************************************************************
430 /* Compression of curve CRVOLD in a table of */
431 /* coeff. of even : CRVNEW(*,0,*) */
432 /* and uneven range : CRVNEW(*,1,*). */
436 /* COMPRESSION,CURVE. */
438 /* INPUT ARGUMENTS : */
439 /* ------------------ */
440 /* NDIM : Space Dimension. */
441 /* NCOFMX : Max nb of coeff. of the curve to compress. */
442 /* NCOEFF : Max nb of coeff. of the compressed curve. */
443 /* CRVOLD : The curve (0:NCOFMX-1,NDIM) to compress. */
445 /* OUTPUT ARGUMENTS : */
446 /* ------------------- */
447 /* CRVNEW : Curve compacted in (0:(NCOEFF-1)/2,0,NDIM) (containing
449 /* even terms) and in (0:(NCOEFF-1)/2,1,NDIM) */
450 /* (containing uneven terms). */
453 /* ---------------- */
455 /* REFERENCES CALLED : */
456 /* ----------------------- */
458 /* DESCRIPTION/NOTES/LIMITATIONS : */
459 /* ----------------------------------- */
460 /* This routine is useful to prepare coefficients of a */
461 /* curve in an orthogonal base (Legendre or Jacobi) before */
462 /* calculating the coefficients in the canonical; base [-1,1] by */
464 /* ***********************************************************************
467 /* Name of the routine */
469 /* Parameter adjustments */
470 crvold_dim1 = *ncofmx;
471 crvold_offset = crvold_dim1;
472 crvold -= crvold_offset;
473 crvnew_dim1 = (*ncoeff - 1) / 2 + 1;
474 crvnew_offset = crvnew_dim1 << 1;
475 crvnew -= crvnew_offset;
478 ibb = AdvApp2Var_SysBase::mnfndeb_();
480 AdvApp2Var_SysBase::mgenmsg_("MMAPCMP", 7L);
483 ndegre = *ncoeff - 1;
485 for (nd = 1; nd <= i__1; ++nd) {
488 for (idg = 0; idg <= i__2; ++idg) {
489 crvnew[idg + (nd << 1) * crvnew_dim1] = crvold[ipair + nd *
498 i__2 = (ndegre - 1) / 2;
499 for (idg = 0; idg <= i__2; ++idg) {
500 crvnew[idg + ((nd << 1) + 1) * crvnew_dim1] = crvold[impair + nd *
511 /* ---------------------------------- The end ---------------------------
515 AdvApp2Var_SysBase::mgsomsg_("MMAPCMP", 7L);
520 //=======================================================================
521 //function : mmaper0_
523 //=======================================================================
524 int mmaper0_(integer *ncofmx,
533 /* System generated locals */
534 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
537 /* Local variables */
542 /* ***********************************************************************
547 /* Calculate the max error of approximation done when */
548 /* only the first NCFNEW coefficients of a curve are preserved.
550 /* Degree NCOEFF-1 written in the base of Legendre (Jacobi */
555 /* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
557 /* INPUT ARGUMENTS : */
558 /* ------------------ */
559 /* NCOFMX : Max. degree of the curve. */
560 /* NDIMEN : Space dimension. */
561 /* NCOEFF : Degree +1 of the curve. */
562 /* CRVLGD : Curve the degree which of should be lowered. */
563 /* NCFNEW : Degree +1 of the resulting polynom. */
565 /* OUTPUT ARGUMENTS : */
566 /* ------------------- */
567 /* YCVMAX : Auxiliary Table (max error on each dimension).
569 /* ERRMAX : Precision of the approximation. */
572 /* ---------------- */
574 /* REFERENCES CALLED : */
575 /* ----------------------- */
577 /* DESCRIPTION/NOTES/LIMITATIONS : */
578 /* ----------------------------------- */
579 /* ***********************************************************************
583 /* ------------------- Init to calculate an error -----------------------
586 /* Parameter adjustments */
588 crvlgd_dim1 = *ncofmx;
589 crvlgd_offset = crvlgd_dim1 + 1;
590 crvlgd -= crvlgd_offset;
594 for (ii = 1; ii <= i__1; ++ii) {
599 /* ------ Minimum that can be reached : Stop at 1 or NCFNEW ------
603 if (*ncfnew + 1 > ncut) {
607 /* -------------- Elimination of high degree coefficients-----------
609 /* ----------- Loop on the series of Legendre: NCUT --> NCOEFF --------
613 for (ii = ncut; ii <= i__1; ++ii) {
614 /* Factor of renormalization (Maximum of Li(t)). */
615 bidon = ((ii - 1) * 2. + 1.) / 2.;
619 for (nd = 1; nd <= i__2; ++nd) {
620 ycvmax[nd] += (d__1 = crvlgd[ii + nd * crvlgd_dim1], advapp_abs(d__1)) *
627 /* -------------- The error is the norm of the vector error ---------------
630 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
632 /* --------------------------------- Fin --------------------------------
638 //=======================================================================
639 //function : mmaper2_
641 //=======================================================================
642 int mmaper2_(integer *ncofmx,
651 /* Initialized data */
653 static doublereal xmaxj[57] = { .9682458365518542212948163499456,
654 .986013297183269340427888048593603,
655 1.07810420343739860362585159028115,
656 1.17325804490920057010925920756025,
657 1.26476561266905634732910520370741,
658 1.35169950227289626684434056681946,
659 1.43424378958284137759129885012494,
660 1.51281316274895465689402798226634,
661 1.5878364329591908800533936587012,
662 1.65970112228228167018443636171226,
663 1.72874345388622461848433443013543,
664 1.7952515611463877544077632304216,
665 1.85947199025328260370244491818047,
666 1.92161634324190018916351663207101,
667 1.98186713586472025397859895825157,
668 2.04038269834980146276967984252188,
669 2.09730119173852573441223706382076,
670 2.15274387655763462685970799663412,
671 2.20681777186342079455059961912859,
672 2.25961782459354604684402726624239,
673 2.31122868752403808176824020121524,
674 2.36172618435386566570998793688131,
675 2.41117852396114589446497298177554,
676 2.45964731268663657873849811095449,
677 2.50718840313973523778244737914028,
678 2.55385260994795361951813645784034,
679 2.59968631659221867834697883938297,
680 2.64473199258285846332860663371298,
681 2.68902863641518586789566216064557,
682 2.73261215675199397407027673053895,
683 2.77551570192374483822124304745691,
684 2.8177699459714315371037628127545,
685 2.85940333797200948896046563785957,
686 2.90044232019793636101516293333324,
687 2.94091151970640874812265419871976,
688 2.98083391718088702956696303389061,
689 3.02023099621926980436221568258656,
690 3.05912287574998661724731962377847,
691 3.09752842783622025614245706196447,
692 3.13546538278134559341444834866301,
693 3.17295042316122606504398054547289,
694 3.2099992681699613513775259670214,
695 3.24662674946606137764916854570219,
696 3.28284687953866689817670991319787,
697 3.31867291347259485044591136879087,
698 3.35411740487202127264475726990106,
699 3.38919225660177218727305224515862,
700 3.42390876691942143189170489271753,
701 3.45827767149820230182596660024454,
702 3.49230918177808483937957161007792,
703 3.5260130200285724149540352829756,
704 3.55939845146044235497103883695448,
705 3.59247431368364585025958062194665,
706 3.62524904377393592090180712976368,
707 3.65773070318071087226169680450936,
708 3.68992700068237648299565823810245,
709 3.72184531357268220291630708234186 };
711 /* System generated locals */
712 integer crvjac_dim1, crvjac_offset, i__1, i__2;
715 /* Local variables */
722 /* ***********************************************************************
727 /* Calculate max approximation error i faite lorsque l' on */
728 /* ne conserve que les premiers NCFNEW coefficients d' une courbe
730 /* de degre NCOEFF-1 ecrite dans la base de Jacobi d' ordre 2. */
734 /* JACOBI, POLYGON, APPROXIMATION, ERROR. */
736 /* INPUT ARGUMENTS : */
737 /* ------------------ */
738 /* NCOFMX : Max. degree of the curve. */
739 /* NDIMEN : Space dimension. */
740 /* NCOEFF : Degree +1 of the curve. */
741 /* CRVLGD : Curve the degree which of should be lowered. */
742 /* NCFNEW : Degree +1 of the resulting polynom. */
744 /* OUTPUT ARGUMENTS : */
745 /* ------------------- */
746 /* YCVMAX : Auxiliary Table (max error on each dimension).
748 /* ERRMAX : Precision of the approximation. */
751 /* ---------------- */
753 /* REFERENCES CALLED : */
754 /* ----------------------- */
755 /* DESCRIPTION/NOTES/LIMITATIONS : */
756 /* ----------------------------------- */
760 /* ------------------ Table of maximums of (1-t2)*Ji(t) ----------------
763 /* Parameter adjustments */
765 crvjac_dim1 = *ncofmx;
766 crvjac_offset = crvjac_dim1 + 1;
767 crvjac -= crvjac_offset;
773 /* ------------------- Init for error calculation -----------------------
777 for (ii = 1; ii <= i__1; ++ii) {
782 /* ------ Min. Degree that can be attained : Stop at 3 or NCFNEW ------
787 i__1 = idec, i__2 = *ncfnew + 1;
788 ncut = advapp_max(i__1,i__2);
790 /* -------------- Removal of coefficients of high degree -----------
792 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
796 for (ii = ncut; ii <= i__1; ++ii) {
797 /* Factor of renormalization. */
798 bidon = xmaxj[ii - idec];
800 for (nd = 1; nd <= i__2; ++nd) {
801 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
808 /* -------------- The error is the norm of the vector error ---------------
811 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
813 /* --------------------------------- Fin --------------------------------
819 /* MAPER4.f -- translated by f2c (version 19960827).
820 You must link the resulting object file with the libraries:
821 -lf2c -lm (in that order)
825 //=======================================================================
826 //function : mmaper4_
828 //=======================================================================
829 int mmaper4_(integer *ncofmx,
837 /* Initialized data */
839 static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
840 1.05299572648705464724876659688996,
841 1.0949715351434178709281698645813,
842 1.15078388379719068145021100764647,
843 1.2094863084718701596278219811869,
844 1.26806623151369531323304177532868,
845 1.32549784426476978866302826176202,
846 1.38142537365039019558329304432581,
847 1.43575531950773585146867625840552,
848 1.48850442653629641402403231015299,
849 1.53973611681876234549146350844736,
850 1.58953193485272191557448229046492,
851 1.63797820416306624705258190017418,
852 1.68515974143594899185621942934906,
853 1.73115699602477936547107755854868,
854 1.77604489805513552087086912113251,
855 1.81989256661534438347398400420601,
856 1.86276344480103110090865609776681,
857 1.90471563564740808542244678597105,
858 1.94580231994751044968731427898046,
859 1.98607219357764450634552790950067,
860 2.02556989246317857340333585562678,
861 2.06433638992049685189059517340452,
862 2.10240936014742726236706004607473,
863 2.13982350649113222745523925190532,
864 2.17661085564771614285379929798896,
865 2.21280102016879766322589373557048,
866 2.2484214321456956597803794333791,
867 2.28349755104077956674135810027654,
868 2.31805304852593774867640120860446,
869 2.35210997297725685169643559615022,
870 2.38568889602346315560143377261814,
871 2.41880904328694215730192284109322,
872 2.45148841120796359750021227795539,
873 2.48374387161372199992570528025315,
874 2.5155912654873773953959098501893,
875 2.54704548720896557684101746505398,
876 2.57812056037881628390134077704127,
877 2.60882970619319538196517982945269,
878 2.63918540521920497868347679257107,
879 2.66919945330942891495458446613851,
880 2.69888301230439621709803756505788,
881 2.72824665609081486737132853370048,
882 2.75730041251405791603760003778285,
883 2.78605380158311346185098508516203,
884 2.81451587035387403267676338931454,
885 2.84269522483114290814009184272637,
886 2.87060005919012917988363332454033,
887 2.89823818258367657739520912946934,
888 2.92561704377132528239806135133273,
889 2.95274375377994262301217318010209,
890 2.97962510678256471794289060402033,
891 3.00626759936182712291041810228171,
892 3.03267744830655121818899164295959,
893 3.05886060707437081434964933864149 };
895 /* System generated locals */
896 integer crvjac_dim1, crvjac_offset, i__1, i__2;
899 /* Local variables */
906 /* ***********************************************************************
911 /* Calculate the max. error of approximation made when */
912 /* only first NCFNEW coefficients of a curve are preserved
914 /* degree NCOEFF-1 is written in the base of Jacobi of order 4. */
917 /* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
919 /* INPUT ARGUMENTS : */
920 /* ------------------ */
921 /* NCOFMX : Max. degree of the curve. */
922 /* NDIMEN : Space dimension. */
923 /* NCOEFF : Degree +1 of the curve. */
924 /* CRVJAC : Curve the degree which of should be lowered. */
925 /* NCFNEW : Degree +1 of the resulting polynom. */
927 /* OUTPUT ARGUMENTS : */
928 /* ------------------- */
929 /* YCVMAX : Auxiliary Table (max error on each dimension).
931 /* ERRMAX : Precision of the approximation. */
934 /* ---------------- */
936 /* REFERENCES CALLED : */
937 /* ----------------------- */
939 /* DESCRIPTION/NOTES/LIMITATIONS : */
942 /* ***********************************************************************
946 /* ---------------- Table of maximums of ((1-t2)2)*Ji(t) ---------------
949 /* Parameter adjustments */
951 crvjac_dim1 = *ncofmx;
952 crvjac_offset = crvjac_dim1 + 1;
953 crvjac -= crvjac_offset;
959 /* ------------------- Init for error calculation -----------------------
963 for (ii = 1; ii <= i__1; ++ii) {
968 /* ------ Min. Degree that can be attained : Stop at 5 or NCFNEW ------
973 i__1 = idec, i__2 = *ncfnew + 1;
974 ncut = advapp_max(i__1,i__2);
976 /* -------------- Removal of high degree coefficients -----------
978 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
982 for (ii = ncut; ii <= i__1; ++ii) {
983 /* Factor of renormalisation. */
984 bidon = xmaxj[ii - idec];
986 for (nd = 1; nd <= i__2; ++nd) {
987 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
994 /* -------------- The error is the norm of the error vector ---------------
997 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
999 /* --------------------------------- End --------------------------------
1005 //=======================================================================
1006 //function : mmaper6_
1008 //=======================================================================
1009 int mmaper6_(integer *ncofmx,
1018 /* Initialized data */
1020 static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
1021 1.11626917091567929907256116528817,
1022 1.1327140810290884106278510474203,
1023 1.1679452722668028753522098022171,
1024 1.20910611986279066645602153641334,
1025 1.25228283758701572089625983127043,
1026 1.29591971597287895911380446311508,
1027 1.3393138157481884258308028584917,
1028 1.3821288728999671920677617491385,
1029 1.42420414683357356104823573391816,
1030 1.46546895108549501306970087318319,
1031 1.50590085198398789708599726315869,
1032 1.54550385142820987194251585145013,
1033 1.58429644271680300005206185490937,
1034 1.62230484071440103826322971668038,
1035 1.65955905239130512405565733793667,
1036 1.69609056468292429853775667485212,
1037 1.73193098017228915881592458573809,
1038 1.7671112206990325429863426635397,
1039 1.80166107681586964987277458875667,
1040 1.83560897003644959204940535551721,
1041 1.86898184653271388435058371983316,
1042 1.90180515174518670797686768515502,
1043 1.93410285411785808749237200054739,
1044 1.96589749778987993293150856865539,
1045 1.99721027139062501070081653790635,
1046 2.02806108474738744005306947877164,
1047 2.05846864831762572089033752595401,
1048 2.08845055210580131460156962214748,
1049 2.11802334209486194329576724042253,
1050 2.14720259305166593214642386780469,
1051 2.17600297710595096918495785742803,
1052 2.20443832785205516555772788192013,
1053 2.2325216999457379530416998244706,
1054 2.2602654243075083168599953074345,
1055 2.28768115912702794202525264301585,
1056 2.3147799369092684021274946755348,
1057 2.34157220782483457076721300512406,
1058 2.36806787963276257263034969490066,
1059 2.39427635443992520016789041085844,
1060 2.42020656255081863955040620243062,
1061 2.44586699364757383088888037359254,
1062 2.47126572552427660024678584642791,
1063 2.49641045058324178349347438430311,
1064 2.52130850028451113942299097584818,
1065 2.54596686772399937214920135190177,
1066 2.5703922285006754089328998222275,
1067 2.59459096001908861492582631591134,
1068 2.61856915936049852435394597597773,
1069 2.64233265984385295286445444361827,
1070 2.66588704638685848486056711408168,
1071 2.68923766976735295746679957665724,
1072 2.71238965987606292679677228666411 };
1074 /* System generated locals */
1075 integer crvjac_dim1, crvjac_offset, i__1, i__2;
1078 /* Local variables */
1085 /* ***********************************************************************
1089 /* Calculate the max. error of approximation made when */
1090 /* only first NCFNEW coefficients of a curve are preserved
1092 /* degree NCOEFF-1 is written in the base of Jacobi of order 6. */
1095 /* JACOBI,POLYGON,APPROXIMATION,ERROR. */
1097 /* INPUT ARGUMENTS : */
1098 /* ------------------ */
1099 /* NCOFMX : Max. degree of the curve. */
1100 /* NDIMEN : Space dimension. */
1101 /* NCOEFF : Degree +1 of the curve. */
1102 /* CRVJAC : Curve the degree which of should be lowered. */
1103 /* NCFNEW : Degree +1 of the resulting polynom. */
1105 /* OUTPUT ARGUMENTS : */
1106 /* ------------------- */
1107 /* YCVMAX : Auxiliary Table (max error on each dimension).
1109 /* ERRMAX : Precision of the approximation. */
1111 /* COMMONS USED : */
1112 /* ---------------- */
1114 /* REFERENCES CALLED : */
1115 /* ----------------------- */
1117 /* DESCRIPTION/NOTES/LIMITATIONS : */
1119 /* ***********************************************************************
1123 /* ---------------- Table of maximums of ((1-t2)3)*Ji(t) ---------------
1126 /* Parameter adjustments */
1128 crvjac_dim1 = *ncofmx;
1129 crvjac_offset = crvjac_dim1 + 1;
1130 crvjac -= crvjac_offset;
1136 /* ------------------- Init for error calculation -----------------------
1140 for (ii = 1; ii <= i__1; ++ii) {
1145 /* ------ Min Degree that can be attained : Stop at 3 or NCFNEW ------
1150 i__1 = idec, i__2 = *ncfnew + 1;
1151 ncut = advapp_max(i__1,i__2);
1153 /* -------------- Removal of high degree coefficients -----------
1155 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
1159 for (ii = ncut; ii <= i__1; ++ii) {
1160 /* Factor of renormalization. */
1161 bidon = xmaxj[ii - idec];
1163 for (nd = 1; nd <= i__2; ++nd) {
1164 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
1171 /* -------------- The error is the norm of the vector error ---------------
1174 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
1176 /* --------------------------------- END --------------------------------
1182 //=======================================================================
1183 //function : AdvApp2Var_MathBase::mmaperx_
1185 //=======================================================================
1186 int AdvApp2Var_MathBase::mmaperx_(integer *ncofmx,
1197 /* System generated locals */
1198 integer crvjac_dim1, crvjac_offset;
1200 /* Local variables */
1203 /* **********************************************************************
1207 /* Calculate the max. error of approximation made when */
1208 /* only first NCFNEW coefficients of a curve are preserved
1210 /* degree NCOEFF-1 is written in the base of Jacobi of order IORDRE. */
1213 /* JACOBI,LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
1215 /* INPUT ARGUMENTS : */
1216 /* ------------------ */
1217 /* NCOFMX : Max. degree of the curve. */
1218 /* NDIMEN : Space dimension. */
1219 /* NCOEFF : Degree +1 of the curve. */
1220 /* IORDRE : Order of continuity at the extremities. */
1221 /* CRVJAC : Curve the degree which of should be lowered. */
1222 /* NCFNEW : Degree +1 of the resulting polynom. */
1224 /* OUTPUT ARGUMENTS : */
1225 /* ------------------- */
1226 /* YCVMAX : Auxiliary Table (max error on each dimension).
1228 /* ERRMAX : Precision of the approximation. */
1229 /* IERCOD = 0, OK */
1230 /* = 1, order of constraints (IORDRE) is not within the */
1231 /* autorized values. */
1232 /* COMMONS USED : */
1233 /* ---------------- */
1235 /* REFERENCES CALLED : */
1236 /* ----------------------- */
1238 /* DESCRIPTION/NOTES/LIMITATIONS : */
1239 /* ----------------------------------- */
1240 /* Canceled and replaced MMAPERR. */
1241 /* ***********************************************************************
1245 /* Parameter adjustments */
1247 crvjac_dim1 = *ncofmx;
1248 crvjac_offset = crvjac_dim1 + 1;
1249 crvjac -= crvjac_offset;
1253 /* --> Order of Jacobi polynoms */
1254 jord = ( *iordre + 1) << 1;
1257 mmaper0_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1259 } else if (jord == 2) {
1260 mmaper2_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1262 } else if (jord == 4) {
1263 mmaper4_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1265 } else if (jord == 6) {
1266 mmaper6_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1272 /* ----------------------------------- Fin ------------------------------
1278 //=======================================================================
1279 //function : mmarc41_
1281 //=======================================================================
1282 int mmarc41_(integer *ndimax,
1292 /* System generated locals */
1293 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
1296 /* Local variables */
1298 doublereal tbaux[61];
1304 /* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
1305 /* IMPLICIT INTEGER (I-N) */
1307 /* ***********************************************************************
1312 /* Creation of curve C2(v) defined on (0,1) identic to */
1313 /* curve C1(u) defined on (U0,U1) (change of parameter */
1318 /* LIMITATION, RESTRICTION, CURVE */
1320 /* INPUT ARGUMENTS : */
1321 /* ------------------ */
1322 /* NDIMAX : Space Dimensioning. */
1323 /* NDIMEN : Curve Dimension. */
1324 /* NCOEFF : Nb of coefficients of the curve. */
1325 /* CRVOLD : Curve to be limited. */
1326 /* UPARA0 : Min limit of the interval limiting the curve.
1328 /* UPARA1 : Max limit of the interval limiting the curve.
1331 /* OUTPUT ARGUMENTS : */
1332 /* ------------------- */
1333 /* CRVNEW : Relimited curve, defined on (0,1) and equal to */
1334 /* CRVOLD defined on (U0,U1). */
1335 /* IERCOD : = 0, OK */
1336 /* =10, Nb of coeff. <1 or > 61. */
1338 /* COMMONS USED : */
1339 /* ---------------- */
1340 /* REFERENCES CALLED : */
1341 /* ---------------------- */
1343 /* MAERMSG MCRFILL MVCVIN2 */
1346 /* DESCRIPTION/NOTES/LIMITATIONS : */
1347 /* ----------------------------------- */
1348 /* ---> Algorithm used in this general case is based on the */
1349 /* following principle : */
1350 /* Let S(t) = a0 + a1*t + a2*t**2 + ... of degree NCOEFF-1, and */
1351 /* U(t) = b0 + b1*t, then the coeff. of */
1352 /* S(U(t)) are calculated step by step with help of table TBAUX. */
1353 /* At each step number N (N=2 to NCOEFF), TBAUX(n) contains */
1354 /* the n-th coefficient of U(t)**N for n=1 to N. (RBD) */
1355 /* ---> Reference : KNUTH, 'The Art of Computer Programming', */
1356 /* Vol. 2/'Seminumerical Algorithms', */
1357 /* Ex. 11 p:451 et solution p:562. (RBD) */
1359 /* ---> Removal of the input argument CRVOLD by CRVNEW is */
1360 /* possible, which means that the call : */
1361 /* CALL MMARC41(NDIMAX,NDIMEN,NCOEFF,CURVE,UPARA0,UPARA1 */
1362 /* ,CURVE,IERCOD) */
1363 /* is absolutely LEGAL. (RBD) */
1366 /* **********************************************************************
1369 /* Name of the routine */
1371 /* Auxiliary table of coefficients of (UPARA1-UPARA0)T+UPARA0 */
1372 /* with power N=1 to NCOEFF-1. */
1375 /* Parameter adjustments */
1376 crvnew_dim1 = *ndimax;
1377 crvnew_offset = crvnew_dim1 + 1;
1378 crvnew -= crvnew_offset;
1379 crvold_dim1 = *ndimax;
1380 crvold_offset = crvold_dim1 + 1;
1381 crvold -= crvold_offset;
1385 /* **********************************************************************
1387 /* CASE WHEN PROCESSING CAN'T BE DONE */
1388 /* **********************************************************************
1390 if (*ncoeff > 61 || *ncoeff < 1) {
1394 /* **********************************************************************
1397 /* **********************************************************************
1399 if (*ndimen == *ndimax && *upara0 == 0. && *upara1 == 1.) {
1400 nboct = (*ndimax << 3) * *ncoeff;
1401 AdvApp2Var_SysBase::mcrfill_(&nboct,
1402 &crvold[crvold_offset],
1403 &crvnew[crvnew_offset]);
1406 /* **********************************************************************
1408 /* INVERSION 3D : FAST PROCESSING */
1409 /* **********************************************************************
1411 if (*upara0 == 1. && *upara1 == 0.) {
1412 if (*ndimen == 3 && *ndimax == 3 && *ncoeff <= 21) {
1413 mvcvinv_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
1417 /* ******************************************************************
1419 /* INVERSION 2D : FAST PROCESSING */
1420 /* ******************************************************************
1422 if (*ndimen == 2 && *ndimax == 2 && *ncoeff <= 21) {
1423 mvcvin2_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
1428 /* **********************************************************************
1430 /* GENERAL PROCESSING */
1431 /* **********************************************************************
1433 /* -------------------------- Initializations ---------------------------
1437 for (nd = 1; nd <= i__1; ++nd) {
1438 crvnew[nd + crvnew_dim1] = crvold[nd + crvold_dim1];
1445 tbaux[1] = *upara1 - *upara0;
1447 /* ----------------------- Calculation of coeff. of CRVNEW ------------------
1451 for (ncf = 2; ncf <= i__1; ++ncf) {
1453 /* ------------ Take into account NCF-th coeff. of CRVOLD --------
1457 for (ncj = 1; ncj <= i__2; ++ncj) {
1458 bid = tbaux[ncj - 1];
1460 for (nd = 1; nd <= i__3; ++nd) {
1461 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + ncf *
1468 bid = tbaux[ncf - 1];
1470 for (nd = 1; nd <= i__2; ++nd) {
1471 crvnew[nd + ncf * crvnew_dim1] = crvold[nd + ncf * crvold_dim1] *
1476 /* --------- Calculate (NCF+1) coeff. of ((U1-U0)*t + U0)**(NCF) ---
1479 bid = *upara1 - *upara0;
1480 tbaux[ncf] = tbaux[ncf - 1] * bid;
1481 for (ncj = ncf; ncj >= 2; --ncj) {
1482 tbaux[ncj - 1] = tbaux[ncj - 1] * *upara0 + tbaux[ncj - 2] * bid;
1485 tbaux[0] *= *upara0;
1490 /* -------------- Take into account the last coeff. of CRVOLD -----------
1494 for (ncj = 1; ncj <= i__1; ++ncj) {
1495 bid = tbaux[ncj - 1];
1497 for (nd = 1; nd <= i__2; ++nd) {
1498 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + *ncoeff *
1505 for (nd = 1; nd <= i__1; ++nd) {
1506 crvnew[nd + *ncoeff * crvnew_dim1] = crvold[nd + *ncoeff *
1507 crvold_dim1] * tbaux[*ncoeff - 1];
1511 /* ---------------------------- The end ---------------------------------
1516 AdvApp2Var_SysBase::maermsg_("MMARC41", iercod, 7L);
1522 //=======================================================================
1523 //function : AdvApp2Var_MathBase::mmarcin_
1525 //=======================================================================
1526 int AdvApp2Var_MathBase::mmarcin_(integer *ndimax,
1536 /* System generated locals */
1537 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
1541 /* Local variables */
1544 doublereal tabaux[61];
1553 /* **********************************************************************
1556 /* Creation of curve C2(v) defined on [U0,U1] identic to */
1557 /* curve C1(u) defined on [-1,1] (change of parameter */
1558 /* of a curve) with INVERSION of indices of the resulting table. */
1562 /* GENERALIZED LIMITATION, RESTRICTION, INVERSION, CURVE */
1564 /* INPUT ARGUMENTS : */
1565 /* ------------------ */
1566 /* NDIMAX : Maximum Space Dimensioning. */
1567 /* NDIMEN : Curve Dimension. */
1568 /* NCOEFF : Nb of coefficients of the curve. */
1569 /* CRVOLD : Curve to be limited. */
1570 /* U0 : Min limit of the interval limiting the curve.
1572 /* U1 : Max limit of the interval limiting the curve.
1575 /* OUTPUT ARGUMENTS : */
1576 /* ------------------- */
1577 /* CRVNEW : Relimited curve, defined on [U0,U1] and equal to */
1578 /* CRVOLD defined on [-1,1]. */
1579 /* IERCOD : = 0, OK */
1580 /* =10, Nb of coeff. <1 or > 61. */
1581 /* =13, the requested interval of variation is null. */
1582 /* COMMONS USED : */
1583 /* ---------------- */
1584 /* REFERENCES CALLED : */
1585 /* ---------------------- */
1586 /* DESCRIPTION/NOTES/LIMITATIONS : */
1587 /* ----------------------------------- */
1589 /* **********************************************************************
1592 /* Name of the routine */
1594 /* Auxiliary table of coefficients of X1*T+X0 */
1595 /* with power N=1 to NCOEFF-1. */
1598 /* Parameter adjustments */
1599 crvnew_dim1 = *ndimax;
1600 crvnew_offset = crvnew_dim1 + 1;
1601 crvnew -= crvnew_offset;
1602 crvold_dim1 = *ncoeff;
1603 crvold_offset = crvold_dim1 + 1;
1604 crvold -= crvold_offset;
1607 ibb = AdvApp2Var_SysBase::mnfndeb_();
1609 AdvApp2Var_SysBase::mgenmsg_("MMARCIN", 7L);
1612 /* At zero machine it is tested if the output interval is not null */
1614 AdvApp2Var_MathBase::mmveps3_(&eps3);
1615 if ((d__1 = *u1 - *u0, advapp_abs(d__1)) < eps3) {
1621 /* **********************************************************************
1623 /* CASE WHEN THE PROCESSING IS IMPOSSIBLE */
1624 /* **********************************************************************
1626 if (*ncoeff > 61 || *ncoeff < 1) {
1630 /* **********************************************************************
1632 /* IF NO CHANGE OF THE INTERVAL OF DEFINITION */
1633 /* (ONLY INVERSION OF INDICES OF TABLE CRVOLD) */
1634 /* **********************************************************************
1636 if (*ndim == *ndimax && *u0 == -1. && *u1 == 1.) {
1637 AdvApp2Var_MathBase::mmcvinv_(ndim, ncoeff, ndim, &crvold[crvold_offset], &crvnew[
1641 /* **********************************************************************
1643 /* CASE WHEN THE NEW INTERVAL OF DEFINITION IS [0,1] */
1644 /* **********************************************************************
1646 if (*u0 == 0. && *u1 == 1.) {
1647 mmcvstd_(ncoeff, ndimax, ncoeff, ndim, &crvold[crvold_offset], &
1648 crvnew[crvnew_offset]);
1651 /* **********************************************************************
1653 /* GENERAL PROCESSING */
1654 /* **********************************************************************
1656 /* -------------------------- Initialization ---------------------------
1659 x0 = -(*u1 + *u0) / (*u1 - *u0);
1660 x1 = 2. / (*u1 - *u0);
1662 for (nd = 1; nd <= i__1; ++nd) {
1663 crvnew[nd + crvnew_dim1] = crvold[nd * crvold_dim1 + 1];
1672 /* ----------------------- Calculation of coeff. of CRVNEW ------------------
1676 for (ncf = 2; ncf <= i__1; ++ncf) {
1678 /* ------------ Take into account the NCF-th coeff. of CRVOLD --------
1682 for (ncj = 1; ncj <= i__2; ++ncj) {
1683 bid = tabaux[ncj - 1];
1685 for (nd = 1; nd <= i__3; ++nd) {
1686 crvnew[nd + ncj * crvnew_dim1] += crvold[ncf + nd *
1693 bid = tabaux[ncf - 1];
1695 for (nd = 1; nd <= i__2; ++nd) {
1696 crvnew[nd + ncf * crvnew_dim1] = crvold[ncf + nd * crvold_dim1] *
1701 /* --------- Calculation of (NCF+1) coeff. of [X1*t + X0]**(NCF) --------
1704 tabaux[ncf] = tabaux[ncf - 1] * x1;
1705 for (ncj = ncf; ncj >= 2; --ncj) {
1706 tabaux[ncj - 1] = tabaux[ncj - 1] * x0 + tabaux[ncj - 2] * x1;
1714 /* -------------- Take into account the last coeff. of CRVOLD -----------
1718 for (ncj = 1; ncj <= i__1; ++ncj) {
1719 bid = tabaux[ncj - 1];
1721 for (nd = 1; nd <= i__2; ++nd) {
1722 crvnew[nd + ncj * crvnew_dim1] += crvold[*ncoeff + nd *
1729 for (nd = 1; nd <= i__1; ++nd) {
1730 crvnew[nd + *ncoeff * crvnew_dim1] = crvold[*ncoeff + nd *
1731 crvold_dim1] * tabaux[*ncoeff - 1];
1735 /* ---------------------------- The end ---------------------------------
1740 AdvApp2Var_SysBase::maermsg_("MMARCIN", iercod, 7L);
1743 AdvApp2Var_SysBase::mgsomsg_("MMARCIN", 7L);
1748 //=======================================================================
1749 //function : mmatvec_
1751 //=======================================================================
1752 int mmatvec_(integer *nligne,
1763 /* System generated locals */
1766 /* Local variables */
1768 integer jmin, jmax, i__, j, k;
1773 /* ***********************************************************************
1778 /* Produce vector matrix in form of profile */
1783 /* RESERVE, MATRIX, PRODUCT, VECTOR, PROFILE */
1785 /* INPUT ARGUMENTS : */
1786 /* -------------------- */
1787 /* NLIGNE : Line number of the matrix of constraints */
1788 /* NCOLON : Number of column of the matrix of constraints */
1789 /* GNSTOC: Number of coefficients in the profile of matrix GMATRI */
1791 /* GPOSIT: Table of positioning of terms of storage */
1792 /* GPOSIT(1,I) contains the number of terms-1 on the line I */
1793 /* in the profile of the matrix. */
1794 /* GPOSIT(2,I) contains the index of storage of diagonal term*/
1796 /* GPOSIT(3,I) contains the index of column of the first term of */
1797 /* profile of line I */
1798 /* GNSTOC: Number of coefficients in the profile of matrix */
1800 /* GMATRI : Matrix of constraints in form of profile */
1801 /* VECIN : Input vector */
1802 /* DEBLIG : Line indexusing which the vector matrix is calculated */
1804 /* OUTPUT ARGUMENTS */
1805 /* --------------------- */
1806 /* VECOUT : VECTOR PRODUCT */
1808 /* IERCOD : ERROR CODE */
1811 /* COMMONS USED : */
1812 /* ------------------ */
1815 /* REFERENCES CALLED : */
1816 /* --------------------- */
1819 /* DESCRIPTION/NOTES/LIMITATIONS : */
1820 /* ----------------------------------- */
1822 /* ***********************************************************************
1825 /* ***********************************************************************
1830 /* ***********************************************************************
1832 /* INITIALISATIONS */
1833 /* ***********************************************************************
1836 /* Parameter adjustments */
1843 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1845 AdvApp2Var_SysBase::mgenmsg_("MMATVEC", 7L);
1849 /* ***********************************************************************
1852 /* ***********************************************************************
1854 AdvApp2Var_SysBase::mvriraz_(nligne,
1857 for (i__ = *deblig; i__ <= i__1; ++i__) {
1859 jmin = gposit[i__ * 3 + 3];
1860 jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
1861 aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
1863 for (j = jmin; j <= i__2; ++j) {
1865 somme += gmatri[k] * vecin[j];
1867 vecout[i__] = somme;
1876 /* ***********************************************************************
1878 /* ERROR PROCESSING */
1879 /* ***********************************************************************
1885 /* ***********************************************************************
1887 /* RETURN CALLING PROGRAM */
1888 /* ***********************************************************************
1893 /* ___ DESALLOCATION, ... */
1895 AdvApp2Var_SysBase::maermsg_("MMATVEC", iercod, 7L);
1897 AdvApp2Var_SysBase::mgsomsg_("MMATVEC", 7L);
1903 //=======================================================================
1904 //function : mmbulld_
1906 //=======================================================================
1907 int AdvApp2Var_MathBase::mmbulld_(integer *nbcoln,
1913 /* System generated locals */
1914 integer dtabtr_dim1, dtabtr_offset, i__1, i__2;
1916 /* Local variables */
1919 integer nite1, nite2, nchan, i1, i2;
1921 /* ***********************************************************************
1926 /* Parsing of columns of a table of integers in increasing order */
1929 /* POINT-ENTRY, PARSING */
1930 /* INPUT ARGUMENTS : */
1931 /* -------------------- */
1932 /* - NBCOLN : Number of columns in the table */
1933 /* - NBLIGN : Number of lines in the table */
1934 /* - DTABTR : Table of integers to be parsed */
1935 /* - NUMCLE : Position of the key on the column */
1937 /* OUTPUT ARGUMENTS : */
1938 /* --------------------- */
1939 /* - DTABTR : Parsed table */
1941 /* COMMONS USED : */
1942 /* ------------------ */
1945 /* REFERENCES CALLED : */
1946 /* --------------------- */
1949 /* DESCRIPTION/NOTES/LIMITATIONS : */
1950 /* ----------------------------------- */
1951 /* Particularly performant if the table is almost parsed */
1952 /* In the opposite case it is better to use MVSHELD */
1953 /* ***********************************************************************
1956 /* Parameter adjustments */
1957 dtabtr_dim1 = *nblign;
1958 dtabtr_offset = dtabtr_dim1 + 1;
1959 dtabtr -= dtabtr_offset;
1962 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1964 AdvApp2Var_SysBase::mgenmsg_("MMBULLD", 7L);
1970 /* ***********************************************************************
1973 /* ***********************************************************************
1976 /* ---->ALGORITHM in N^2 / 2 additional iteration */
1980 /* ----> Parsing from left to the right */
1984 for (i1 = nite2; i1 <= i__1; ++i1) {
1985 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1 - 1)
1988 for (i2 = 1; i2 <= i__2; ++i2) {
1989 daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
1990 dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
1992 dtabtr[i2 + i1 * dtabtr_dim1] = daux;
2001 /* ----> Parsing from right to the left */
2006 for (i1 = nite1; i1 >= i__1; --i1) {
2007 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1
2008 - 1) * dtabtr_dim1]) {
2010 for (i2 = 1; i2 <= i__2; ++i2) {
2011 daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
2012 dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
2014 dtabtr[i2 + i1 * dtabtr_dim1] = daux;
2028 /* ***********************************************************************
2030 /* ERROR PROCESSING */
2031 /* ***********************************************************************
2034 /* ----> No errors at calling functions, only tests and loops. */
2036 /* ***********************************************************************
2038 /* RETURN CALLING PROGRAM */
2039 /* ***********************************************************************
2045 AdvApp2Var_SysBase::mgsomsg_("MMBULLD", 7L);
2052 //=======================================================================
2053 //function : AdvApp2Var_MathBase::mmcdriv_
2055 //=======================================================================
2056 int AdvApp2Var_MathBase::mmcdriv_(integer *ndimen,
2065 /* System generated locals */
2066 integer courbe_dim1, courbe_offset, crvdrv_dim1, crvdrv_offset, i__1,
2069 /* Local variables */
2071 doublereal mfactk, bid;
2074 /* ***********************************************************************
2079 /* Calculate matrix of a derivate curve of order IDERIV. */
2080 /* with input parameters other than output parameters. */
2085 /* COEFFICIENTS,CURVE,DERIVATE I-EME. */
2087 /* INPUT ARGUMENTS : */
2088 /* ------------------ */
2089 /* NDIMEN : Space dimension (2 or 3 in general) */
2090 /* NCOEFF : Degree +1 of the curve. */
2091 /* COURBE : Table of coefficients of the curve. */
2092 /* IDERIV : Required order of derivation : 1=1st derivate, etc... */
2094 /* OUTPUT ARGUMENTS : */
2095 /* ------------------- */
2096 /* NCOFDV : Degree +1 of the derivative of order IDERIV of the curve. */
2097 /* CRVDRV : Table of coefficients of the derivative of order IDERIV */
2100 /* COMMONS USED : */
2101 /* ---------------- */
2103 /* REFERENCES CALLED : */
2104 /* ----------------------- */
2106 /* DESCRIPTION/NOTES/LIMITATIONS : */
2107 /* ----------------------------------- */
2109 /* ---> It is possible to take as output argument the curve */
2110 /* and the number of coeff passed at input by making : */
2111 /* CALL MMCDRIV(NDIMEN,NCOEFF,COURBE,IDERIV,NCOEFF,COURBE). */
2112 /* After this call, NCOEFF does the number of coeff of the derived */
2113 /* curve the coefficients which of are stored in CURVE. */
2114 /* Attention to the coefficients of CURVE of rank superior to */
2115 /* NCOEFF : they are not set to zero. */
2117 /* ---> Algorithm : */
2118 /* The code below was written basing on the following algorithm:
2121 /* Let P(t) = a1 + a2*t + ... an*t**n. Derivate of order k of P */
2122 /* (containing n-k coefficients) is calculated as follows : */
2124 /* Pk(t) = a(k+1)*CNP(k,k)*k! */
2125 /* + a(k+2)*CNP(k+1,k)*k! * t */
2129 /* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
2130 /* ***********************************************************************
2134 /* -------------- Case when the order of derivative is -------------------
2136 /* ---------------- greater than the degree of the curve ---------------------
2139 /* **********************************************************************
2144 /* Serves to provide the coefficients of binome (Pascal's triangle). */
2148 /* Binomial coeff from 0 to 60. read only . init par block data */
2150 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
2151 /* ----------------------------------- */
2152 /* Binomial coefficients form a triangular matrix. */
2153 /* This matrix is completed in table CNP by its transposition. */
2154 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
2156 /* Initialization is done by block-data MMLLL09.RES, */
2157 /* created by program MQINICNP.FOR). */
2158 /* **********************************************************************
2163 /* ***********************************************************************
2166 /* Parameter adjustments */
2167 crvdrv_dim1 = *ndimen;
2168 crvdrv_offset = crvdrv_dim1 + 1;
2169 crvdrv -= crvdrv_offset;
2170 courbe_dim1 = *ndimen;
2171 courbe_offset = courbe_dim1 + 1;
2172 courbe -= courbe_offset;
2175 if (*ideriv >= *ncoeff) {
2177 for (i__ = 1; i__ <= i__1; ++i__) {
2178 crvdrv[i__ + crvdrv_dim1] = 0.;
2184 /* **********************************************************************
2186 /* General processing */
2187 /* **********************************************************************
2189 /* --------------------- Calculation of Factorial(IDERIV) ------------------
2195 for (i__ = 2; i__ <= i__1; ++i__) {
2200 /* ------------ Calculation of coeff of the derived of order IDERIV ----------
2202 /* ---> Attention : coefficient binomial C(n,m) is represented in */
2203 /* MCCNP by CNP(N+1,M+1). */
2206 for (j = k + 1; j <= i__1; ++j) {
2207 bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
2209 for (i__ = 1; i__ <= i__2; ++i__) {
2210 crvdrv[i__ + (j - k) * crvdrv_dim1] = bid * courbe[i__ + j *
2217 *ncofdv = *ncoeff - *ideriv;
2219 /* -------------------------------- The end -----------------------------
2226 //=======================================================================
2227 //function : AdvApp2Var_MathBase::mmcglc1_
2229 //=======================================================================
2230 int AdvApp2Var_MathBase::mmcglc1_(integer *ndimax,
2243 /* System generated locals */
2244 integer courbe_dim1, courbe_offset, i__1;
2247 /* Local variables */
2249 doublereal tdeb, tfin;
2251 doublereal oldso = 0.;
2255 doublereal dif, pas;
2259 /* ***********************************************************************
2264 /* Allows calculating the length of an arc of curve POLYNOMIAL */
2265 /* on an interval [A,B]. */
2269 /* LENGTH,CURVE,GAUSS,PRIVATE. */
2271 /* INPUT ARGUMENTS : */
2272 /* ------------------ */
2273 /* NDIMAX : Max. number of lines of tables */
2274 /* (i.e. max. nb of polynoms). */
2275 /* NDIMEN : Dimension of the space (nb of polynoms). */
2276 /* NCOEFF : Nb of coefficients of the polynom. This is degree + 1.
2278 /* COURBE(NDIMAX,NCOEFF) : Coefficients of the curve. */
2279 /* TDEBUT : Lower limit of the interval of integration for */
2280 /* length calculation. */
2281 /* TFINAL : Upper limit of the interval of integration for */
2282 /* length calculation. */
2283 /* EPSILN : REQIRED precision for length calculation. */
2285 /* OUTPUT ARGUMENTS : */
2286 /* ------------------- */
2287 /* XLONGC : Length of the arc of curve */
2288 /* ERREUR : Precision OBTAINED for the length calculation. */
2289 /* IERCOD : Error code, 0 OK, >0 Serious error. */
2290 /* = 1 Too much iterations, the best calculated resultat */
2291 /* (is almost ERROR) */
2292 /* = 2 Pb MMLONCV (no result) */
2293 /* = 3 NDIM or NCOEFF invalid (no result) */
2295 /* COMMONS USED : */
2296 /* ---------------- */
2298 /* REFERENCES CALLED : */
2299 /* ----------------------- */
2301 /* DESCRIPTION/NOTES/LIMITATIONS : */
2302 /* ----------------------------------- */
2303 /* The polynom is actually a set of polynoms with */
2304 /* coefficients arranged in a table of 2 indices, */
2305 /* each line relative to the polynom. */
2306 /* The polynom is defined by these coefficients ordered */
2307 /* by increasing power of the variable. */
2308 /* All polynoms have the same number of coefficients (the */
2311 /* This program cancels and replaces LENGCV, MLONGC and MLENCV. */
2313 /* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
2316 /* ***********************************************************************
2319 /* Name of the routine */
2322 /* ------------------------ General Initialization ---------------------
2325 /* Parameter adjustments */
2326 courbe_dim1 = *ndimax;
2327 courbe_offset = courbe_dim1 + 1;
2328 courbe -= courbe_offset;
2331 ibb = AdvApp2Var_SysBase::mnfndeb_();
2333 AdvApp2Var_SysBase::mgenmsg_("MMCGLC1", 7L);
2340 /* ------ Test of equity of limits */
2342 if (*tdebut == *tfinal) {
2347 /* ------ Test of the dimension and the number of coefficients */
2349 if (*ndimen <= 0 || *ncoeff <= 0) {
2353 /* ----- Nb of current cutting, nb of iteration, */
2354 /* max nb of iterations */
2361 /* ------ Variation of the nb of intervals */
2362 /* Multiplied by 2 at each iteration */
2365 pas = (*tfinal - *tdebut) / ndec;
2368 /* ------ Loop on all current NDEC intervals */
2371 for (kk = 1; kk <= i__1; ++kk) {
2373 /* ------ Limits of the current integration interval */
2375 tdeb = *tdebut + (kk - 1) * pas;
2377 mmloncv_(ndimax, ndimen, ncoeff, &courbe[courbe_offset], &tdeb, &tfin,
2389 /* ----------------- Test of the maximum number of iterations ------------
2392 /* Test if passes at least once ** */
2401 /* ------ Take into account DIF - Test of convergence */
2404 dif = (d__1 = sottc - oldso, advapp_abs(d__1));
2406 /* ------ If DIF is OK, leave..., otherwise: */
2408 if (dif > *epsiln) {
2410 /* ------ If nb iteration exceeded, leave */
2417 /* ------ Otherwise continue by cutting the initial interval.
2427 /* ------------------------------ THE END -------------------------------
2435 /* ---> PB in MMLONCV */
2441 /* ---> NCOEFF or NDIM invalid. */
2449 AdvApp2Var_SysBase::maermsg_("MMCGLC1", iercod, 7L);
2452 AdvApp2Var_SysBase::mgsomsg_("MMCGLC1", 7L);
2457 //=======================================================================
2458 //function : mmchole_
2460 //=======================================================================
2461 int mmchole_(integer *,//mxcoef,
2470 /* System generated locals */
2471 integer i__1, i__2, i__3;
2474 /* Builtin functions */
2477 /* Local variables */
2479 integer kmin, i__, j, k;
2481 integer ptini, ptcou;
2484 /* ***********************************************************************
2489 /* Produce decomposition of choleski of matrix A in S.S */
2490 /* Calculate inferior triangular matrix S. */
2494 /* RESOLUTION, MFACTORISATION, MATRIX_PROFILE, CHOLESKI */
2496 /* INPUT ARGUMENTS : */
2497 /* -------------------- */
2498 /* MXCOEF : Max number of terms in the hessian profile */
2499 /* DIMENS : Dimension of the problem */
2500 /* AMATRI(MXCOEF) : Coefficients of the matrix profile */
2501 /* APOSIT(1,*) : Distance diagonal-left extremity of the line
2503 /* APOSIT(2,*) : Position of diagonal terms in HESSIE */
2504 /* POSUIV(MXCOEF) : first line inferior not out of profile */
2506 /* OUTPUT ARGUMENTS : */
2507 /* --------------------- */
2508 /* CHOMAT(MXCOEF) : Inferior triangular matrix preserving the */
2509 /* profile of AMATRI. */
2510 /* IERCOD : error code */
2512 /* = 1 : non-defined positive matrix */
2514 /* COMMONS USED : */
2515 /* ------------------ */
2519 /* REFERENCES CALLED : */
2520 /* ---------------------- */
2522 /* DESCRIPTION/NOTES/LIMITATIONS : */
2523 /* ----------------------------------- */
2524 /* DEBUG LEVEL = 4 */
2525 /* ***********************************************************************
2528 /* ***********************************************************************
2533 /* ***********************************************************************
2535 /* INITIALISATIONS */
2536 /* ***********************************************************************
2539 /* Parameter adjustments */
2546 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
2548 AdvApp2Var_SysBase::mgenmsg_("MMCHOLE", 7L);
2552 /* ***********************************************************************
2555 /* ***********************************************************************
2559 for (j = 1; j <= i__1; ++j) {
2561 ptini = aposit[(j << 1) + 2];
2565 for (k = ptini - aposit[(j << 1) + 1]; k <= i__2; ++k) {
2566 /* Computing 2nd power */
2568 somme += d__1 * d__1;
2571 if (amatri[ptini] - somme < 1e-32) {
2574 chomat[ptini] = sqrt(amatri[ptini] - somme);
2578 while(posuiv[ptcou] > 0) {
2580 i__ = posuiv[ptcou];
2581 ptcou = aposit[(i__ << 1) + 2] - (i__ - j);
2583 /* Calculate the sum of S .S for k =1 a j-1 */
2587 i__2 = i__ - aposit[(i__ << 1) + 1], i__3 = j - aposit[(j << 1) +
2589 kmin = advapp_max(i__2,i__3);
2591 for (k = kmin; k <= i__2; ++k) {
2592 somme += chomat[aposit[(i__ << 1) + 2] - (i__ - k)] * chomat[
2593 aposit[(j << 1) + 2] - (j - k)];
2596 chomat[ptcou] = (amatri[ptcou] - somme) / chomat[ptini];
2602 /* ***********************************************************************
2604 /* ERROR PROCESSING */
2605 /* ***********************************************************************
2612 /* ***********************************************************************
2614 /* RETURN CALLING PROGRAM */
2615 /* ***********************************************************************
2620 AdvApp2Var_SysBase::maermsg_("MMCHOLE", iercod, 7L);
2622 AdvApp2Var_SysBase::mgsomsg_("MMCHOLE", 7L);
2628 //=======================================================================
2629 //function : AdvApp2Var_MathBase::mmcvctx_
2631 //=======================================================================
2632 int AdvApp2Var_MathBase::mmcvctx_(integer *ndimen,
2642 /* System generated locals */
2643 integer ctrtes_dim1, ctrtes_offset, crvres_dim1, crvres_offset,
2644 xmatri_dim1, xmatri_offset, tabaux_dim1, tabaux_offset, i__1,
2647 /* Local variables */
2648 integer moup1, nordr;
2650 integer ibb, ncf, ndv;
2654 /* ***********************************************************************
2659 /* Calculate a polynomial curve checking the */
2660 /* passage constraints (interpolation) */
2661 /* from first derivatives, etc... to extremities. */
2662 /* Parameters at the extremities are supposed to be -1 and 1. */
2666 /* ALL, AB_SPECIFI::CONSTRAINTS&,INTERPOLATION,&CURVE */
2668 /* INPUT ARGUMENTS : */
2669 /* ------------------ */
2670 /* NDIMEN : Space Dimension. */
2671 /* NCOFMX : Nb of coeff. of curve CRVRES on each */
2673 /* NDERIV : Order of constraint with derivatives : */
2674 /* 0 --> interpolation simple. */
2675 /* 1 --> interpolation+constraints with 1st. */
2676 /* 2 --> cas (0)+ (1) + " " 2nd derivatives. */
2678 /* CTRTES : Table of constraints. */
2679 /* CTRTES(*,1,*) = contraints at -1. */
2680 /* CTRTES(*,2,*) = contraints at 1. */
2682 /* OUTPUT ARGUMENTS : */
2683 /* ------------------- */
2684 /* CRVRES : Resulting curve defined on (-1,1). */
2685 /* TABAUX : Auxilliary matrix. */
2686 /* XMATRI : Auxilliary matrix. */
2688 /* COMMONS UTILISES : */
2689 /* ---------------- */
2693 /* REFERENCES CALLED : */
2694 /* ---------------------- */
2696 /* MAERMSG R*8 DFLOAT MGENMSG */
2697 /* MGSOMSG MMEPS1 MMRSLW */
2700 /* DESCRIPTION/NOTES/LIMITATIONS : */
2701 /* ----------------------------------- */
2702 /* The polynom (or the curve) is calculated by solving a */
2703 /* system of linear equations. If the imposed degree is great */
2704 /* it is preferable to call a routine based on */
2705 /* Lagrange or Hermite interpolation depending on the case. */
2706 /* (for a high degree the matrix of the system can be badly */
2707 /* conditionned). */
2708 /* This routine returns a curve defined in (-1,1). */
2709 /* In general case, it is necessary to use MCVCTG. */
2711 /* ***********************************************************************
2714 /* Name of the routine */
2717 /* Parameter adjustments */
2718 crvres_dim1 = *ncofmx;
2719 crvres_offset = crvres_dim1 + 1;
2720 crvres -= crvres_offset;
2721 xmatri_dim1 = *nderiv + 1;
2722 xmatri_offset = xmatri_dim1 + 1;
2723 xmatri -= xmatri_offset;
2724 tabaux_dim1 = *nderiv + 1 + *ndimen;
2725 tabaux_offset = tabaux_dim1 + 1;
2726 tabaux -= tabaux_offset;
2727 ctrtes_dim1 = *ndimen;
2728 ctrtes_offset = ctrtes_dim1 * 3 + 1;
2729 ctrtes -= ctrtes_offset;
2732 ibb = AdvApp2Var_SysBase::mnfndeb_();
2734 AdvApp2Var_SysBase::mgenmsg_("MMCVCTX", 7L);
2737 AdvApp2Var_MathBase::mmeps1_(&eps1);
2739 /* ****************** CALCULATION OF EVEN COEFFICIENTS *********************
2741 /* ------------------------- Initialization -----------------------------
2744 nordr = *nderiv + 1;
2746 for (ncf = 1; ncf <= i__1; ++ncf) {
2747 tabaux[ncf + tabaux_dim1] = 1.;
2751 /* ---------------- Calculation of terms corresponding to derivatives -------
2755 for (ndv = 2; ndv <= i__1; ++ndv) {
2757 for (ncf = 1; ncf <= i__2; ++ncf) {
2758 tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
2759 tabaux_dim1] * (doublereal) ((ncf << 1) - ndv);
2765 /* ------------------ Writing the second member -----------------------
2770 for (ndv = 1; ndv <= i__1; ++ndv) {
2772 for (nd = 1; nd <= i__2; ++nd) {
2773 tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
2774 + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2775 * ctrtes_dim1]) / 2.;
2782 /* -------------------- Resolution of the system ---------------------------
2785 mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2786 xmatri_offset], iercod);
2791 for (nd = 1; nd <= i__1; ++nd) {
2793 for (ncf = 1; ncf <= i__2; ++ncf) {
2794 crvres[(ncf << 1) - 1 + nd * crvres_dim1] = xmatri[ncf + nd *
2801 /* ***************** CALCULATION OF UNEVEN COEFFICIENTS ********************
2803 /* ------------------------- Initialization -----------------------------
2808 for (ncf = 1; ncf <= i__1; ++ncf) {
2809 tabaux[ncf + tabaux_dim1] = 1.;
2813 /* ---------------- Calculation of terms corresponding to derivatives -------
2817 for (ndv = 2; ndv <= i__1; ++ndv) {
2819 for (ncf = 1; ncf <= i__2; ++ncf) {
2820 tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
2821 tabaux_dim1] * (doublereal) ((ncf << 1) - ndv + 1);
2827 /* ------------------ Writing of the second member -----------------------
2832 for (ndv = 1; ndv <= i__1; ++ndv) {
2834 for (nd = 1; nd <= i__2; ++nd) {
2835 tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
2836 + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2837 * ctrtes_dim1]) / 2.;
2844 /* -------------------- Solution of the system ---------------------------
2847 mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2848 xmatri_offset], iercod);
2853 for (nd = 1; nd <= i__1; ++nd) {
2855 for (ncf = 1; ncf <= i__2; ++ncf) {
2856 crvres[(ncf << 1) + nd * crvres_dim1] = xmatri[ncf + nd *
2863 /* --------------------------- The end ----------------------------------
2868 AdvApp2Var_SysBase::maermsg_("MMCVCTX", iercod, 7L);
2871 AdvApp2Var_SysBase::mgsomsg_("MMCVCTX", 7L);
2877 //=======================================================================
2878 //function : AdvApp2Var_MathBase::mmcvinv_
2880 //=======================================================================
2881 int AdvApp2Var_MathBase::mmcvinv_(integer *ndimax,
2888 /* Initialized data */
2890 static char nomprg[8+1] = "MMCVINV ";
2892 /* System generated locals */
2893 integer curve_dim1, curve_offset, curveo_dim1, curveo_offset, i__1, i__2;
2895 /* Local variables */
2896 integer i__, nd, ibb;
2899 /* ***********************************************************************
2904 /* Inversion of arguments of the final curve. */
2908 /* SMOOTHING,CURVE */
2911 /* INPUT ARGUMENTS : */
2912 /* ------------------ */
2914 /* NDIM: Space Dimension. */
2915 /* NCOEF: Degree of the polynom. */
2916 /* CURVEO: The curve before inversion. */
2918 /* OUTPUT ARGUMENTS : */
2919 /* ------------------- */
2920 /* CURVE: The curve after inversion. */
2922 /* COMMONS USED : */
2923 /* ---------------- */
2924 /* REFERENCES APPELEES : */
2925 /* ----------------------- */
2926 /* DESCRIPTION/NOTES/LIMITATIONS : */
2927 /* ----------------------------------- */
2928 /* ***********************************************************************
2931 /* The name of the routine */
2932 /* Parameter adjustments */
2933 curve_dim1 = *ndimax;
2934 curve_offset = curve_dim1 + 1;
2935 curve -= curve_offset;
2936 curveo_dim1 = *ncoef;
2937 curveo_offset = curveo_dim1 + 1;
2938 curveo -= curveo_offset;
2942 ibb = AdvApp2Var_SysBase::mnfndeb_();
2944 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
2948 for (i__ = 1; i__ <= i__1; ++i__) {
2950 for (nd = 1; nd <= i__2; ++nd) {
2951 curve[nd + i__ * curve_dim1] = curveo[i__ + nd * curveo_dim1];
2960 //=======================================================================
2961 //function : mmcvstd_
2963 //=======================================================================
2964 int mmcvstd_(integer *ncofmx,
2972 /* System generated locals */
2973 integer courbe_dim1, crvcan_dim1, crvcan_offset, i__1, i__2, i__3;
2975 /* Local variables */
2976 integer ndeg, i__, j, j1, nd, ibb;
2980 /* ***********************************************************************
2985 /* Transform curve defined between [-1,1] into [0,1]. */
2989 /* LIMITATION,RESTRICTION,CURVE */
2991 /* INPUT ARGUMENTS : */
2992 /* ------------------ */
2993 /* NDIMAX : Dimension of the space. */
2994 /* NDIMEN : Dimension of the curve. */
2995 /* NCOEFF : Degree of the curve. */
2996 /* CRVCAN(NCOFMX,NDIMEN): The curve is defined at the interval [-1,1]. */
2998 /* OUTPUT ARGUMENTS : */
2999 /* ------------------- */
3000 /* CURVE(NDIMAX,NCOEFF): Curve defined at the interval [0,1]. */
3002 /* COMMONS USED : */
3003 /* ---------------- */
3005 /* REFERENCES CALLED : */
3006 /* ----------------------- */
3008 /* DESCRIPTION/NOTES/LIMITATIONS : */
3009 /* ----------------------------------- */
3011 /* ***********************************************************************
3014 /* Name of the program. */
3017 /* **********************************************************************
3022 /* Provides binomial coefficients (Pascal triangle). */
3026 /* Binomial coefficient from 0 to 60. read only . init by block data */
3028 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3029 /* ----------------------------------- */
3030 /* Binomial coefficients form a triangular matrix. */
3031 /* This matrix is completed in table CNP by its transposition. */
3032 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
3034 /* Initialization is done with block-data MMLLL09.RES, */
3035 /* created by the program MQINICNP.FOR. */
3037 /* **********************************************************************
3042 /* ***********************************************************************
3045 /* Parameter adjustments */
3046 courbe_dim1 = *ndimax;
3048 crvcan_dim1 = *ncofmx;
3049 crvcan_offset = crvcan_dim1;
3050 crvcan -= crvcan_offset;
3053 ibb = AdvApp2Var_SysBase::mnfndeb_();
3055 AdvApp2Var_SysBase::mgenmsg_("MMCVSTD", 7L);
3059 /* ------------------ Construction of the resulting curve ----------------
3063 for (nd = 1; nd <= i__1; ++nd) {
3065 for (j = 0; j <= i__2; ++j) {
3068 for (i__ = j; i__ <= i__3; i__ += 2) {
3069 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
3073 courbe[nd + j * courbe_dim1] = bid;
3078 for (i__ = j1; i__ <= i__3; i__ += 2) {
3079 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
3083 courbe[nd + j * courbe_dim1] -= bid;
3089 /* ------------------- Renormalization of the CURVE -------------------------
3094 for (i__ = 0; i__ <= i__1; ++i__) {
3096 for (nd = 1; nd <= i__2; ++nd) {
3097 courbe[nd + i__ * courbe_dim1] *= bid;
3104 /* ----------------------------- The end --------------------------------
3108 AdvApp2Var_SysBase::mgsomsg_("MMCVSTD", 7L);
3113 //=======================================================================
3114 //function : AdvApp2Var_MathBase::mmdrc11_
3116 //=======================================================================
3117 int AdvApp2Var_MathBase::mmdrc11_(integer *iordre,
3122 doublereal *mfactab)
3125 /* System generated locals */
3126 integer courbe_dim1, courbe_offset, points_dim2, points_offset, i__1,
3129 /* Local variables */
3131 integer ndeg, i__, j, ndgcb, nd, ibb;
3134 /* **********************************************************************
3139 /* Calculation of successive derivatives of equation CURVE with */
3140 /* parameters -1, 1 from order 0 to order IORDRE */
3141 /* included. The calculation is produced without knowing the coefficients of */
3142 /* derivatives of the curve. */
3146 /* POSITIONING,EXTREMITIES,CURVE,DERIVATIVE. */
3148 /* INPUT ARGUMENTS : */
3149 /* ------------------ */
3150 /* IORDRE : Maximum order of calculation of derivatives. */
3151 /* NDIMEN : Dimension of the space. */
3152 /* NCOEFF : Number of coefficients of the curve (degree+1). */
3153 /* COURBE : Table of coefficients of the curve. */
3155 /* OUTPUT ARGUMENTS : */
3156 /* ------------------- */
3157 /* POINTS : Table of values of consecutive derivatives */
3158 /* of parameters -1.D0 and 1.D0. */
3159 /* MFACTAB : Auxiliary table for calculation of factorial(I).
3162 /* COMMONS USED : */
3163 /* ---------------- */
3166 /* REFERENCES CALLED : */
3167 /* ----------------------- */
3169 /* DESCRIPTION/NOTES/LIMITATIONS : */
3170 /* ----------------------------------- */
3172 /* ---> ATTENTION, the coefficients of the curve are */
3173 /* in a reverse order. */
3175 /* ---> The algorithm of calculation of derivatives is based on */
3176 /* generalization of Horner scheme : */
3178 /* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
3181 /* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
3183 /* aj = a(j-1).x + u(k-j) */
3184 /* bj = b(j-1).x + a(j-1) */
3185 /* cj = c(j-1).x + b(j-1) */
3187 /* So : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
3189 /* The algorithm is generalized easily for calculation of */
3196 /* Reference : D. KNUTH, "The Art of Computer Programming" */
3197 /* --------- Vol. 2/Seminumerical Algorithms */
3198 /* Addison-Wesley Pub. Co. (1969) */
3199 /* pages 423-425. */
3201 /* **********************************************************************
3204 /* Name of the routine */
3206 /* Parameter adjustments */
3207 points_dim2 = *iordre + 1;
3208 points_offset = (points_dim2 << 1) + 1;
3209 points -= points_offset;
3210 courbe_dim1 = *ncoeff;
3211 courbe_offset = courbe_dim1;
3212 courbe -= courbe_offset;
3215 ibb = AdvApp2Var_SysBase::mnfndeb_();
3217 AdvApp2Var_SysBase::mgenmsg_("MMDRC11", 7L);
3220 if (*iordre < 0 || *ncoeff < 1) {
3224 /* ------------------- Initialization of table POINTS -----------------
3227 ndgcb = *ncoeff - 1;
3229 for (nd = 1; nd <= i__1; ++nd) {
3230 points[(nd * points_dim2 << 1) + 1] = courbe[ndgcb + nd * courbe_dim1]
3232 points[(nd * points_dim2 << 1) + 2] = courbe[ndgcb + nd * courbe_dim1]
3238 for (nd = 1; nd <= i__1; ++nd) {
3240 for (j = 1; j <= i__2; ++j) {
3241 points[((j + nd * points_dim2) << 1) + 1] = 0.;
3242 points[((j + nd * points_dim2) << 1) + 2] = 0.;
3248 /* Calculation with parameter -1 and 1 */
3251 for (nd = 1; nd <= i__1; ++nd) {
3253 for (ndeg = 1; ndeg <= i__2; ++ndeg) {
3254 for (i__ = *iordre; i__ >= 1; --i__) {
3255 points[((i__ + nd * points_dim2) << 1) + 1] = -points[((i__ + nd
3256 * points_dim2) << 1) + 1] + points[((i__ - 1 + nd *
3257 points_dim2) << 1) + 1];
3258 points[((i__ + nd * points_dim2) << 1) + 2] += points[((i__ - 1
3259 + nd * points_dim2) << 1) + 2];
3262 points[(nd * points_dim2 << 1) + 1] = -points[(nd * points_dim2 <<
3263 1) + 1] + courbe[ndgcb - ndeg + nd * courbe_dim1];
3264 points[(nd * points_dim2 << 1) + 2] += courbe[ndgcb - ndeg + nd *
3271 /* --------------------- Multiplication by factorial(I) --------------
3275 mfac_(&mfactab[1], iordre);
3278 for (nd = 1; nd <= i__1; ++nd) {
3280 for (i__ = 2; i__ <= i__2; ++i__) {
3281 points[((i__ + nd * points_dim2) << 1) + 1] = mfactab[i__] *
3282 points[((i__ + nd * points_dim2) << 1) + 1];
3283 points[((i__ + nd * points_dim2) << 1) + 2] = mfactab[i__] *
3284 points[((i__ + nd * points_dim2) << 1) + 2];
3291 /* ---------------------------- End -------------------------------------
3296 AdvApp2Var_SysBase::mgsomsg_("MMDRC11", 7L);
3301 //=======================================================================
3302 //function : mmdrvcb_
3304 //=======================================================================
3305 int mmdrvcb_(integer *ideriv,
3314 /* System generated locals */
3315 integer courbe_dim1, tabpnt_dim1, i__1, i__2, i__3;
3317 /* Local variables */
3318 integer ndeg, i__, j, nd, ndgcrb, iptpnt, ibb;
3321 /* *********************************************************************** */
3325 /* Calculation of successive derivatives of equation CURVE with */
3326 /* parameter TPARAM from order 0 to order IDERIV included. */
3327 /* The calculation is produced without knowing the coefficients of */
3328 /* derivatives of the CURVE. */
3332 /* POSITIONING,PARAMETER,CURVE,DERIVATIVE. */
3334 /* INPUT ARGUMENTS : */
3335 /* ------------------ */
3336 /* IORDRE : Maximum order of calculation of derivatives. */
3337 /* NDIMEN : Dimension of the space. */
3338 /* NCOEFF : Number of coefficients of the curve (degree+1). */
3339 /* COURBE : Table of coefficients of the curve. */
3340 /* TPARAM : Value of the parameter where the curve should be evaluated. */
3342 /* OUTPUT ARGUMENTS : */
3343 /* ------------------- */
3344 /* TABPNT : Table of values of consecutive derivatives */
3345 /* of parameter TPARAM. */
3346 /* IERCOD : 0 = OK, */
3347 /* 1 = incoherent input. */
3349 /* COMMONS USED : */
3350 /* ---------------- */
3353 /* REFERENCES CALLED : */
3354 /* ----------------------- */
3356 /* DESCRIPTION/NOTES/LIMITATIONS : */
3357 /* ----------------------------------- */
3359 /* The algorithm of calculation of derivatives is based on */
3360 /* generalization of the Horner scheme : */
3362 /* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
3365 /* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
3367 /* aj = a(j-1).x + u(k-j) */
3368 /* bj = b(j-1).x + a(j-1) */
3369 /* cj = c(j-1).x + b(j-1) */
3371 /* So, it is obtained : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
3373 /* The algorithm can be easily generalized for the calculation of */
3380 /* Reference : D. KNUTH, "The Art of Computer Programming" */
3381 /* --------- Vol. 2/Seminumerical Algorithms */
3382 /* Addison-Wesley Pub. Co. (1969) */
3383 /* pages 423-425. */
3385 /* ---> To evaluare derivatives at 0 and 1, it is preferable */
3386 /* to use routine MDRV01.FOR . */
3388 /* **********************************************************************
3391 /* Name of the routine */
3393 /* Parameter adjustments */
3394 tabpnt_dim1 = *ndim;
3396 courbe_dim1 = *ndim;
3400 ibb = AdvApp2Var_SysBase::mnfndeb_();
3402 AdvApp2Var_SysBase::mgenmsg_("MMDRVCB", 7L);
3405 if (*ideriv < 0 || *ncoeff < 1) {
3411 /* ------------------- Initialization of table TABPNT -----------------
3414 ndgcrb = *ncoeff - 1;
3416 for (nd = 1; nd <= i__1; ++nd) {
3417 tabpnt[nd] = courbe[nd + ndgcrb * courbe_dim1];
3424 iptpnt = *ndim * *ideriv;
3425 AdvApp2Var_SysBase::mvriraz_(&iptpnt,
3426 &tabpnt[tabpnt_dim1 + 1]);
3429 /* ------------------------ Calculation of parameter TPARAM ------------------
3433 for (ndeg = 1; ndeg <= i__1; ++ndeg) {
3435 for (nd = 1; nd <= i__2; ++nd) {
3436 for (i__ = *ideriv; i__ >= 1; --i__) {
3437 tabpnt[nd + i__ * tabpnt_dim1] = tabpnt[nd + i__ *
3438 tabpnt_dim1] * *tparam + tabpnt[nd + (i__ - 1) *
3442 tabpnt[nd] = tabpnt[nd] * *tparam + courbe[nd + (ndgcrb - ndeg) *
3449 /* --------------------- Multiplication by factorial(I) -------------
3453 for (i__ = 2; i__ <= i__1; ++i__) {
3455 for (j = 2; j <= i__2; ++j) {
3457 for (nd = 1; nd <= i__3; ++nd) {
3458 tabpnt[nd + i__ * tabpnt_dim1] = (doublereal) j * tabpnt[nd +
3467 /* --------------------------- The end ---------------------------------
3472 AdvApp2Var_SysBase::maermsg_("MMDRVCB", iercod, 7L);
3477 //=======================================================================
3478 //function : AdvApp2Var_MathBase::mmdrvck_
3480 //=======================================================================
3481 int AdvApp2Var_MathBase::mmdrvck_(integer *ncoeff,
3489 /* Initialized data */
3491 static doublereal mmfack[21] = { 1.,2.,6.,24.,120.,720.,5040.,40320.,
3492 362880.,3628800.,39916800.,479001600.,6227020800.,87178291200.,
3493 1.307674368e12,2.0922789888e13,3.55687428096e14,6.402373705728e15,
3494 1.21645100408832e17,2.43290200817664e18,5.109094217170944e19 };
3496 /* System generated locals */
3497 integer courbe_dim1, courbe_offset, i__1, i__2;
3499 /* Local variables */
3500 integer i__, j, k, nd;
3501 doublereal mfactk, bid;
3504 /* IMPLICIT INTEGER (I-N) */
3505 /* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
3508 /* ***********************************************************************
3513 /* Calculate the value of a derived curve of order IDERIV in */
3514 /* a point of parameter TPARAM. */
3518 /* POSITIONING,CURVE,DERIVATIVE of ORDER K. */
3520 /* INPUT ARGUMENTS : */
3521 /* ------------------ */
3522 /* NCOEFF : Degree +1 of the curve. */
3523 /* NDIMEN : Dimension of the space (2 or 3 in general) */
3524 /* COURBE : Table of coefficients of the curve. */
3525 /* IDERIV : Required order of derivation : 1=1st derivative, etc... */
3526 /* TPARAM : Value of parameter of the curve. */
3528 /* OUTPUT ARGUMENTS : */
3529 /* ------------------- */
3530 /* PNTCRB : Point of parameter TPARAM on the derivative of order */
3531 /* IDERIV of CURVE. */
3533 /* COMMONS USED : */
3534 /* ---------------- */
3537 /* REFERENCES CALLED : */
3538 /* ---------------------- */
3540 /* DESCRIPTION/NOTES/LIMITATIONS : */
3541 /* ----------------------------------- */
3543 /* The code below was written basing on the following algorithm :
3546 /* Let P(t) = a1 + a2*t + ... an*t**n. The derivative of order k of P */
3547 /* (containing n-k coefficients) is calculated as follows : */
3549 /* Pk(t) = a(k+1)*CNP(k,k)*k! */
3550 /* + a(k+2)*CNP(k+1,k)*k! * t */
3554 /* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
3556 /* Evaluation is produced following the classic Horner scheme. */
3558 /* ***********************************************************************
3562 /* Factorials (1 to 21) caculated on VAX in R*16 */
3565 /* **********************************************************************
3570 /* Serves to provide binomial coefficients (Pascal triangle). */
3574 /* Binomial Coeff from 0 to 60. read only . init by block data */
3576 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3577 /* ----------------------------------- */
3578 /* Binomial coefficients form a triangular matrix. */
3579 /* This matrix is completed in table CNP by its transposition. */
3580 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
3582 /* Initialization is done by block-data MMLLL09.RES, */
3583 /* created by program MQINICNP.FOR. */
3585 /* **********************************************************************
3590 /* ***********************************************************************
3593 /* Parameter adjustments */
3595 courbe_dim1 = *ndimen;
3596 courbe_offset = courbe_dim1 + 1;
3597 courbe -= courbe_offset;
3601 /* -------------- Case when the order of derivative is greater than -------------------
3603 /* ---------------- the degree of the curve ---------------------
3606 if (*ideriv >= *ncoeff) {
3608 for (nd = 1; nd <= i__1; ++nd) {
3614 /* **********************************************************************
3616 /* General processing*/
3617 /* **********************************************************************
3619 /* --------------------- Calculation of Factorial(IDERIV) ------------------
3623 if (*ideriv <= 21 && *ideriv > 0) {
3624 mfactk = mmfack[k - 1];
3628 for (i__ = 2; i__ <= i__1; ++i__) {
3634 /* ------- Calculation of derivative of order IDERIV of CURVE in TPARAM -----
3636 /* ---> Attention : binomial coefficient C(n,m) is represented in */
3637 /* MCCNP by CNP(N,M). */
3640 for (nd = 1; nd <= i__1; ++nd) {
3641 pntcrb[nd] = courbe[nd + *ncoeff * courbe_dim1] * mmcmcnp_.cnp[*
3642 ncoeff - 1 + k * 61] * mfactk;
3647 for (j = *ncoeff - 1; j >= i__1; --j) {
3648 bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
3650 for (nd = 1; nd <= i__2; ++nd) {
3651 pntcrb[nd] = pntcrb[nd] * *tparam + courbe[nd + j * courbe_dim1] *
3658 /* -------------------------------- The end -----------------------------
3666 //=======================================================================
3667 //function : AdvApp2Var_MathBase::mmeps1_
3669 //=======================================================================
3670 int AdvApp2Var_MathBase::mmeps1_(doublereal *epsilo)
3673 /* ***********************************************************************
3678 /* Extraction of EPS1 from COMMON MPRCSN. EPS1 is spatial zero */
3679 /* equal to 1.D-9 */
3683 /* MPRCSN,PRECISON,EPS1. */
3685 /* INPUT ARGUMENTS : */
3686 /* ------------------ */
3689 /* OUTPUT ARGUMENTS : */
3690 /* ------------------- */
3691 /* EPSILO : Value of EPS1 (spatial zero (10**-9)) */
3693 /* COMMONS USED : */
3694 /* ---------------- */
3696 /* REFERENCES CALLED : */
3697 /* ----------------------- */
3699 /* DESCRIPTION/NOTES/LIMITATIONS : */
3700 /* ----------------------------------- */
3701 /* EPS1 is ABSOLUTE spatial zero, so it is necessary */
3702 /* to use it whenever it is necessary to test if a variable */
3703 /* is null. For example, if the norm of a vector is lower than */
3704 /* EPS1, this vector is NULL ! (when one works in */
3705 /* REAL*8) It is absolutely not advised to test arguments */
3706 /* compared to EPS1**2. Taking into account the rounding errors inevitable */
3707 /* during calculations, this causes testing compared to 0.D0. */
3709 /* ***********************************************************************
3714 /* ***********************************************************************
3719 /* Gives tolerances of invalidity in stream */
3720 /* as well as limits of iterative processes */
3722 /* general context, modifiable by the user */
3726 /* PARAMETER , TOLERANCE */
3728 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3729 /* ----------------------------------- */
3730 /* INITIALISATION : profile , **VIA MPRFTX** at input in stream */
3731 /* loading of default values of the profile in MPRFTX at input */
3732 /* in stream. They are preserved in local variables of MPRFTX */
3734 /* Reset of default values : MDFINT */
3735 /* Interactive modification by the user : MDBINT */
3737 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
3738 /* MEPSPB ... EPS3,EPS4 */
3739 /* MEPSLN ... EPS2, NITERM , NITERR */
3740 /* MEPSNR ... EPS2 , NITERM */
3741 /* MITERR ... NITERR */
3743 /* ***********************************************************************
3746 /* NITERM : max nb of iterations */
3747 /* NITERR : nb of rapid iterations */
3748 /* EPS1 : tolerance of 3D null distance */
3749 /* EPS2 : tolerance of parametric null distance */
3750 /* EPS3 : tolerance to avoid division by 0.. */
3751 /* EPS4 : angular tolerance */
3755 /* ***********************************************************************
3757 *epsilo = mmprcsn_.eps1;
3762 //=======================================================================
3763 //function : mmexthi_
3765 //=======================================================================
3766 int mmexthi_(integer *ndegre,
3770 /* System generated locals */
3773 /* Local variables */
3774 integer iadd, ideb, ndeg2, nmod2, ii, ibb;
3777 /* **********************************************************************
3782 /* Extract of common LDGRTL the weight of formulas of */
3783 /* Gauss quadrature on all roots of Legendre polynoms of degree */
3784 /* NDEGRE defined on [-1,1]. */
3788 /* ALL, AB_SPECIFI::COMMON&, EXTRACTION, &WEIGHT, &GAUSS. */
3790 /* INPUT ARGUMENTS : */
3791 /* ------------------ */
3792 /* NDEGRE : Mathematic degree of Legendre polynom. It should have */
3793 /* 2 <= NDEGRE <= 61. */
3795 /* OUTPUT ARGUMENTS : */
3796 /* ------------------- */
3797 /* HWGAUS : The table of weights of Gauss quadrature formulas */
3798 /* relative to NDEGRE roots of a polynome de Legendre de */
3801 /* COMMONS UTILISES : */
3802 /* ---------------- */
3805 /* REFERENCES CALLED : */
3806 /* ----------------------- */
3808 /* DESCRIPTION/NOTES/LIMITATIONS : */
3809 /* ----------------------------------- */
3810 /* ATTENTION: The condition on NDEGRE ( 2 <= NDEGRE <= 61) is not */
3811 /* tested. The caller should make the test. */
3813 /* Name of the routine */
3816 /* Common MLGDRTL: */
3817 /* This common includes POSITIVE roots of Legendre polynims */
3818 /* AND weights of Gauss quadrature formulas on all */
3819 /* POSITIVE roots of Legendre polynoms. */
3823 /* ***********************************************************************
3828 /* The common of Legendre roots. */
3834 /* DESCRIPTION/NOTES/LIMITATIONS : */
3835 /* ----------------------------------- */
3837 /* ***********************************************************************
3843 /* ROOTAB : Table of all roots of Legendre polynoms */
3844 /* within the interval [0,1]. They are ranked for the degrees increasing from */
3846 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
3847 /* The adressing is the same. */
3848 /* HI0TAB : Table of Legendre interpolators for root x=0 */
3849 /* of polynoms of UNEVEN degree. */
3850 /* RTLTB0 : Table of Li(uk) where uk are the roots of */
3851 /* Legendre polynom of EVEN degree. */
3852 /* RTLTB1 : Table of Li(uk) where uk are the roots of */
3853 /* Legendre polynom of UNEVEN degree. */
3856 /************************************************************************
3858 /* Parameter adjustments */
3862 ibb = AdvApp2Var_SysBase::mnfndeb_();
3864 AdvApp2Var_SysBase::mgenmsg_("MMEXTHI", 7L);
3867 ndeg2 = *ndegre / 2;
3868 nmod2 = *ndegre % 2;
3870 /* Address of Gauss weight associated to the 1st strictly */
3871 /* positive root of Legendre polynom of degree NDEGRE in MLGDRTL. */
3873 iadd = ndeg2 * (ndeg2 - 1) / 2 + 1;
3875 /* Index of the 1st HWGAUS element associated to the 1st strictly */
3876 /* positive root of Legendre polynom of degree NDEGRE. */
3878 ideb = (*ndegre + 1) / 2 + 1;
3880 /* Reading of weights associated to strictly positive roots. */
3883 for (ii = ideb; ii <= i__1; ++ii) {
3884 kpt = iadd + ii - ideb;
3885 hwgaus[ii] = mlgdrtl_.hiltab[kpt + nmod2 * 465 - 1];
3889 /* For strictly negative roots, the weight is the same. */
3890 /* i.e HW(1) = HW(NDEGRE), HW(2) = HW(NDEGRE-1), etc... */
3893 for (ii = 1; ii <= i__1; ++ii) {
3894 hwgaus[ii] = hwgaus[*ndegre + 1 - ii];
3898 /* Case of uneven NDEGRE, 0 is root of Legendre polynom, */
3899 /* associated Gauss weights are loaded. */
3902 hwgaus[ndeg2 + 1] = mlgdrtl_.hi0tab[ndeg2];
3905 /* --------------------------- The end ----------------------------------
3909 AdvApp2Var_SysBase::mgsomsg_("MMEXTHI", 7L);
3914 //=======================================================================
3915 //function : mmextrl_
3917 //=======================================================================
3918 int mmextrl_(integer *ndegre,
3921 /* System generated locals */
3924 /* Local variables */
3925 integer iadd, ideb, ndeg2, nmod2, ii, ibb;
3929 /* **********************************************************************
3934 /* Extract of the Common LDGRTL of Legendre polynom roots */
3935 /* of degree NDEGRE defined on [-1,1]. */
3939 /* ALL, AB_SPECIFI::COMMON&, EXTRACTION, &ROOT, &LEGENDRE. */
3941 /* INPUT ARGUMENTS : */
3942 /* ------------------ */
3943 /* NDEGRE : Mathematic degree of Legendre polynom. */
3944 /* It is required to have 2 <= NDEGRE <= 61. */
3946 /* OUTPUT ARGUMENTS : */
3947 /* ------------------- */
3948 /* ROOTLG : The table of roots of Legendre polynom of degree */
3949 /* NDEGRE defined on [-1,1]. */
3951 /* COMMONS USED : */
3952 /* ---------------- */
3955 /* REFERENCES CALLED : */
3956 /* ----------------------- */
3958 /* DESCRIPTION/NOTES/LIMITATIONS : */
3959 /* ----------------------------------- */
3960 /* ATTENTION: Condition of NDEGRE ( 2 <= NDEGRE <= 61) is not */
3961 /* tested. The caller should make the test. */
3963 /* **********************************************************************
3967 /* Name of the routine */
3970 /* Common MLGDRTL: */
3971 /* This common includes POSITIVE roots of Legendre polynoms */
3972 /* AND the weight of Gauss quadrature formulas on all */
3973 /* POSITIVE roots of Legendre polynoms. */
3975 /* ***********************************************************************
3980 /* The common of Legendre roots. */
3987 /* ***********************************************************************
3990 /* ROOTAB : Table of all roots of Legendre polynoms */
3991 /* within the interval [0,1]. They are ranked for the degrees increasing from */
3993 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
3994 /* The adressing is the same. */
3995 /* HI0TAB : Table of Legendre interpolators for root x=0 */
3996 /* of polynoms of UNEVEN degree. */
3997 /* RTLTB0 : Table of Li(uk) where uk are the roots of */
3998 /* Legendre polynom of EVEN degree. */
3999 /* RTLTB1 : Table of Li(uk) where uk are the roots of */
4000 /* Legendre polynom of UNEVEN degree. */
4003 /************************************************************************
4005 /* Parameter adjustments */
4009 ibb = AdvApp2Var_SysBase::mnfndeb_();
4011 AdvApp2Var_SysBase::mgenmsg_("MMEXTRL", 7L);
4014 ndeg2 = *ndegre / 2;
4015 nmod2 = *ndegre % 2;
4017 /* Address of the 1st strictly positive root of Legendre polynom */
4018 /* of degree NDEGRE in MLGDRTL. */
4020 iadd = ndeg2 * (ndeg2 - 1) / 2 + 1;
4022 /* Indice, in ROOTLG, of the 1st strictly positive root */
4023 /* of Legendre polynom of degree NDEGRE. */
4025 ideb = (*ndegre + 1) / 2 + 1;
4027 /* Reading of strictly positive roots. */
4030 for (ii = ideb; ii <= i__1; ++ii) {
4031 kpt = iadd + ii - ideb;
4032 rootlg[ii] = mlgdrtl_.rootab[kpt + nmod2 * 465 - 1];
4036 /* Strictly negative roots are equal to positive roots
4038 /* to the sign i.e RT(1) = -RT(NDEGRE), RT(2) = -RT(NDEGRE-1), etc...
4042 for (ii = 1; ii <= i__1; ++ii) {
4043 rootlg[ii] = -rootlg[*ndegre + 1 - ii];
4047 /* Case NDEGRE uneven, 0 is root of Legendre polynom. */
4050 rootlg[ndeg2 + 1] = 0.;
4053 /* -------------------------------- THE END -----------------------------
4057 AdvApp2Var_SysBase::mgenmsg_("MMEXTRL", 7L);
4062 //=======================================================================
4063 //function : AdvApp2Var_MathBase::mmfmca8_
4065 //=======================================================================
4066 int AdvApp2Var_MathBase::mmfmca8_(const integer *ndimen,
4067 const integer *ncoefu,
4068 const integer *ncoefv,
4069 const integer *ndimax,
4070 const integer *ncfumx,
4071 const integer *,//ncfvmx,
4076 /* System generated locals */
4077 integer tabini_dim1, tabini_dim2, tabini_offset, tabres_dim1, tabres_dim2,
4080 /* Local variables */
4081 integer i__, j, k, ilong;
4085 /* **********************************************************************
4090 /* Expansion of a table containing only most important things into a */
4091 /* greater data table. */
4095 /* ALL, MATH_ACCES:: CARREAU&, DECOMPRESSION, &CARREAU */
4097 /* INPUT ARGUMENTS : */
4098 /* ------------------ */
4099 /* NDIMEN: Dimension of the workspace. */
4100 /* NCOEFU: Degree +1 of the table by u. */
4101 /* NCOEFV: Degree +1 of the table by v. */
4102 /* NDIMAX: Max dimension of the space. */
4103 /* NCFUMX: Max Degree +1 of the table by u. */
4104 /* NCFVMX: Max Degree +1 of the table by v. */
4105 /* TABINI: The table to be decompressed. */
4107 /* OUTPUT ARGUMENTS : */
4108 /* ------------------- */
4109 /* TABRES: Decompressed table. */
4111 /* COMMONS USED : */
4112 /* ---------------- */
4114 /* REFERENCES CALLED : */
4115 /* ----------------------- */
4117 /* DESCRIPTION/NOTES/LIMITATIONS : */
4118 /* ----------------------------------- */
4119 /* The following call : */
4121 /* CALL MMFMCA8(NDIMEN,NCOEFU,NCOEFV,NDIMAX,NCFUMX,NCFVMX,TABINI,TABINI)
4124 /* where TABINI is input/output argument, is possible provided */
4125 /* that the caller has declared TABINI in (NDIMAX,NCFUMX,NCFVMX) */
4127 /* ATTENTION : it is not checked that NDIMAX >= NDIMEN, */
4128 /* NCOEFU >= NCFMXU and NCOEFV >= NCFMXV. */
4130 /* **********************************************************************
4134 /* Parameter adjustments */
4135 tabini_dim1 = *ndimen;
4136 tabini_dim2 = *ncoefu;
4137 tabini_offset = tabini_dim1 * (tabini_dim2 + 1) + 1;
4138 tabini -= tabini_offset;
4139 tabres_dim1 = *ndimax;
4140 tabres_dim2 = *ncfumx;
4141 tabres_offset = tabres_dim1 * (tabres_dim2 + 1) + 1;
4142 tabres -= tabres_offset;
4145 if (*ndimax == *ndimen) {
4149 /* ----------------------- decompression NDIMAX<>NDIMEN -----------------
4152 for (k = *ncoefv; k >= 1; --k) {
4153 for (j = *ncoefu; j >= 1; --j) {
4154 for (i__ = *ndimen; i__ >= 1; --i__) {
4155 tabres[i__ + (j + k * tabres_dim2) * tabres_dim1] = tabini[
4156 i__ + (j + k * tabini_dim2) * tabini_dim1];
4165 /* ----------------------- decompression NDIMAX=NDIMEN ------------------
4169 if (*ncoefu == *ncfumx) {
4172 ilong = (*ndimen << 3) * *ncoefu;
4173 for (k = *ncoefv; k >= 1; --k) {
4174 AdvApp2Var_SysBase::mcrfill_(&ilong,
4175 &tabini[(k * tabini_dim2 + 1) * tabini_dim1 + 1],
4176 &tabres[(k * tabres_dim2 + 1) * tabres_dim1 + 1]);
4181 /* ----------------- decompression NDIMAX=NDIMEN,NCOEFU=NCFUMX ----------
4185 ilong = (*ndimen << 3) * *ncoefu * *ncoefv;
4186 AdvApp2Var_SysBase::mcrfill_(&ilong,
4187 &tabini[tabini_offset],
4188 &tabres[tabres_offset]);
4191 /* ---------------------------- The end ---------------------------------
4198 //=======================================================================
4199 //function : AdvApp2Var_MathBase::mmfmca9_
4201 //=======================================================================
4202 int AdvApp2Var_MathBase::mmfmca9_(integer *ndimax,
4212 /* System generated locals */
4213 integer tabini_dim1, tabini_dim2, tabini_offset, tabres_dim1, tabres_dim2,
4214 tabres_offset, i__1, i__2, i__3;
4216 /* Local variables */
4217 integer i__, j, k, ilong;
4221 /* **********************************************************************
4226 /* Compression of a data table in a table */
4227 /* containing only the main data (the input table is not removed). */
4231 /* ALL, MATH_ACCES:: CARREAU&, COMPRESSION, &CARREAU */
4233 /* INPUT ARGUMENTS : */
4234 /* ------------------ */
4235 /* NDIMAX: Max dimension of the space. */
4236 /* NCFUMX: Max degree +1 of the table by u. */
4237 /* NCFVMX: Max degree +1 of the table by v. */
4238 /* NDIMEN: Dimension of the workspace. */
4239 /* NCOEFU: Degree +1 of the table by u. */
4240 /* NCOEFV: Degree +1 of the table by v. */
4241 /* TABINI: The table to compress. */
4243 /* OUTPUT ARGUMENTS : */
4244 /* ------------------- */
4245 /* TABRES: The compressed table. */
4247 /* COMMONS USED : */
4248 /* ---------------- */
4250 /* REFERENCES CALLED : */
4251 /* ----------------------- */
4253 /* DESCRIPTION/NOTES/LIMITATIONS : */
4254 /* ----------------------------------- */
4255 /* The following call : */
4257 /* CALL MMFMCA9(NDIMAX,NCFUMX,NCFVMX,NDIMEN,NCOEFU,NCOEFV,TABINI,TABINI)
4260 /* where TABINI is input/output argument, is possible provided */
4261 /* that the caller has checked that : */
4263 /* NDIMAX > NDIMEN, */
4264 /* or NDIMAX = NDIMEN and NCFUMX > NCOEFU */
4265 /* or NDIMAX = NDIMEN, NCFUMX = NCOEFU and NCFVMX > NCOEFV */
4267 /* These conditions are not tested in the program. */
4270 /* **********************************************************************
4274 /* Parameter adjustments */
4275 tabini_dim1 = *ndimax;
4276 tabini_dim2 = *ncfumx;
4277 tabini_offset = tabini_dim1 * (tabini_dim2 + 1) + 1;
4278 tabini -= tabini_offset;
4279 tabres_dim1 = *ndimen;
4280 tabres_dim2 = *ncoefu;
4281 tabres_offset = tabres_dim1 * (tabres_dim2 + 1) + 1;
4282 tabres -= tabres_offset;
4285 if (*ndimen == *ndimax) {
4289 /* ----------------------- Compression NDIMEN<>NDIMAX -------------------
4293 for (k = 1; k <= i__1; ++k) {
4295 for (j = 1; j <= i__2; ++j) {
4297 for (i__ = 1; i__ <= i__3; ++i__) {
4298 tabres[i__ + (j + k * tabres_dim2) * tabres_dim1] = tabini[
4299 i__ + (j + k * tabini_dim2) * tabini_dim1];
4308 /* ----------------------- Compression NDIMEN=NDIMAX --------------------
4312 if (*ncoefu == *ncfumx) {
4315 ilong = (*ndimen << 3) * *ncoefu;
4317 for (k = 1; k <= i__1; ++k) {
4318 AdvApp2Var_SysBase::mcrfill_(&ilong,
4319 &tabini[(k * tabini_dim2 + 1) * tabini_dim1 + 1],
4320 &tabres[(k * tabres_dim2 + 1) * tabres_dim1 + 1]);
4325 /* ----------------- Compression NDIMEN=NDIMAX,NCOEFU=NCFUMX ------------
4329 ilong = (*ndimen << 3) * *ncoefu * *ncoefv;
4330 AdvApp2Var_SysBase::mcrfill_(&ilong,
4331 &tabini[tabini_offset],
4332 &tabres[tabres_offset]);
4335 /* ---------------------------- The end ---------------------------------
4342 //=======================================================================
4343 //function : AdvApp2Var_MathBase::mmfmcar_
4345 //=======================================================================
4346 int AdvApp2Var_MathBase::mmfmcar_(integer *ndimen,
4360 /* System generated locals */
4361 integer patold_dim1, patold_dim2, patnew_dim1, patnew_dim2,
4362 i__1, patold_offset,patnew_offset;
4364 /* Local variables */
4365 doublereal* tbaux = 0;
4366 integer ksize, numax, kk;
4370 /* ***********************************************************************
4375 /* LIMITATION OF A SQUARE DEFINED ON (0,1)*(0,1) BETWEEN ISOS */
4376 /* UPARA1 AND UPARA2 (BY U) AND VPARA1 AND VPARA2 BY V. */
4380 /* LIMITATION , SQUARE , PARAMETER */
4382 /* INPUT ARGUMENTS : */
4383 /* ------------------ */
4384 /* NCOFMX: MAX NUMBER OF COEFF OF THE SQUARE BY U */
4385 /* NCOEFU: NUMBER OF COEFF OF THE SQUARE BY U */
4386 /* NCOEFV: NUMBER OF COEFF OF THE SQUARE BY V */
4387 /* PATOLD : THE SQUARE IS LIMITED BY UPARA1,UPARA2 AND VPARA1,VPARA2
4389 /* UPARA1 : LOWER LIMIT OF U */
4390 /* UPARA2 : UPPER LIMIT OF U */
4391 /* VPARA1 : LOWER LIMIT OF V */
4392 /* VPARA2 : UPPER LIMIT OF V */
4394 /* OUTPUT ARGUMENTS : */
4395 /* ------------------- */
4396 /* PATNEW : RELIMITED SQUARE, DEFINED ON (0,1)**2 */
4397 /* IERCOD : =10 COEFF NB TOO GREAT OR NULL */
4398 /* =13 PB IN THE DYNAMIC ALLOCATION */
4401 /* COMMONS USED : */
4402 /* ---------------- */
4404 /* DESCRIPTION/NOTES/LIMITATIONS : */
4405 /* ----------------------------------- */
4406 /* ---> The following call : */
4407 /* CALL MMFMCAR(NCOFMX,NCOEFU,NCOEFV,PATOLD,UPARA1,UPARA2,VPARA1,VPARA2
4410 /* where PATOLD is input/output argument is absolutely legal. */
4412 /* ---> The max number of coeff by u and v of PATOLD is 61 */
4414 /* ---> If NCOEFU < NCOFMX, the data is compressed by MMFMCA9 before */
4415 /* limitation by v to get time during the execution */
4416 /* of MMARC41 that follows (the square is processed as a curve of
4418 /* dimension NDIMEN*NCOEFU possessing NCOEFV coefficients). */
4420 /* ***********************************************************************
4423 /* Name of the routine */
4426 /* Parameter adjustments */
4427 patnew_dim1 = *ndimen;
4428 patnew_dim2 = *ncofmx;
4429 patnew_offset = patnew_dim1 * (patnew_dim2 + 1) + 1;
4430 patnew -= patnew_offset;
4431 patold_dim1 = *ndimen;
4432 patold_dim2 = *ncofmx;
4433 patold_offset = patold_dim1 * (patold_dim2 + 1) + 1;
4434 patold -= patold_offset;
4437 ibb = AdvApp2Var_SysBase::mnfndeb_();
4439 AdvApp2Var_SysBase::mgenmsg_("MMFMCAR", 7L);
4443 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
4445 /* **********************************************************************
4447 /* TEST OF COEFFICIENT NUMBERS */
4448 /* **********************************************************************
4451 if (*ncofmx < *ncoefu) {
4455 if (*ncoefu < 1 || *ncoefu > 61 || *ncoefv < 1 || *ncoefv > 61) {
4460 /* **********************************************************************
4462 /* CASE WHEN UPARA1=VPARA1=0 AND UPARA2=VPARA2=1 */
4463 /* **********************************************************************
4466 if (*upara1 == 0. && *upara2 == 1. && *vpara1 == 0. && *vpara2 == 1.) {
4467 ksize = (*ndimen << 3) * *ncofmx * *ncoefv;
4468 AdvApp2Var_SysBase::mcrfill_(&ksize,
4469 &patold[patold_offset],
4470 &patnew[patnew_offset]);
4474 /* **********************************************************************
4476 /* LIMITATION BY U */
4477 /* **********************************************************************
4480 if (*upara1 == 0. && *upara2 == 1.) {
4484 for (kk = 1; kk <= i__1; ++kk) {
4485 mmarc41_(ndimen, ndimen, ncoefu, &patold[(kk * patold_dim2 + 1) *
4486 patold_dim1 + 1], upara1, upara2, &patnew[(kk * patnew_dim2 +
4487 1) * patnew_dim1 + 1], iercod);
4491 /* **********************************************************************
4493 /* LIMITATION BY V */
4494 /* **********************************************************************
4498 if (*vpara1 == 0. && *vpara2 == 1.) {
4502 /* ----------- LIMITATION BY V (WITH COMPRESSION I.E. NCOEFU<NCOFMX) ----
4505 numax = *ndimen * *ncoefu;
4506 if (*ncofmx != *ncoefu) {
4507 /* ------------------------- Dynamic allocation -------------------
4509 ksize = *ndimen * *ncoefu * *ncoefv;
4510 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &ksize, tbaux, &iofst, &ier);
4515 /* --------------- Compression by (NDIMEN,NCOEFU,NCOEFV) ------------
4517 if (*upara1 == 0. && *upara2 == 1.) {
4518 AdvApp2Var_MathBase::mmfmca9_(ndimen,
4524 &patold[patold_offset],
4527 AdvApp2Var_MathBase::mmfmca9_(ndimen,
4533 &patnew[patnew_offset],
4536 /* ------------------------- Limitation by v ------------------------
4538 mmarc41_(&numax, &numax, ncoefv, &tbaux[iofst], vpara1, vpara2, &
4539 tbaux[iofst], iercod);
4540 /* --------------------- Expansion of TBAUX into PATNEW -------------
4542 AdvApp2Var_MathBase::mmfmca8_(ndimen, ncoefu, ncoefv, ndimen, ncofmx, ncoefv, &tbaux[iofst]
4543 , &patnew[patnew_offset]);
4546 /* -------- LIMITATION BY V (WITHOUT COMPRESSION I.E. NCOEFU=NCOFMX) ---
4550 if (*upara1 == 0. && *upara2 == 1.) {
4551 mmarc41_(&numax, &numax, ncoefv, &patold[patold_offset], vpara1,
4552 vpara2, &patnew[patnew_offset], iercod);
4554 mmarc41_(&numax, &numax, ncoefv, &patnew[patnew_offset], vpara1,
4555 vpara2, &patnew[patnew_offset], iercod);
4560 /* **********************************************************************
4563 /* **********************************************************************
4568 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &ksize, tbaux, &iofst, &ier);
4574 /* ------------------------------ The end -------------------------------
4579 AdvApp2Var_SysBase::maermsg_("MMFMCAR", iercod, 7L);
4582 AdvApp2Var_SysBase::mgsomsg_("MMFMCAR", 7L);
4588 //=======================================================================
4589 //function : AdvApp2Var_MathBase::mmfmcb5_
4591 //=======================================================================
4592 int AdvApp2Var_MathBase::mmfmcb5_(integer *isenmsc,
4603 /* System generated locals */
4604 integer courb1_dim1, courb1_offset, courb2_dim1, courb2_offset, i__1,
4607 /* Local variables */
4608 integer i__, nboct, nd;
4611 /* **********************************************************************
4616 /* Reformating (and eventual compression/decompression) of curve */
4617 /* (ndim,.) by (.,ndim) and vice versa. */
4621 /* ALL , MATH_ACCES :: */
4622 /* COURBE&, REORGANISATION,COMPRESSION,INVERSION , &COURBE */
4624 /* INPUT ARGUMENTS : */
4625 /* -------------------- */
4626 /* ISENMSC : required direction of the transfer : */
4627 /* 1 : passage of (NDIMEN,.) ---> (.,NDIMEN) direction to AB
4629 /* -1 : passage of (.,NDIMEN) ---> (NDIMEN,.) direction to TS,T
4631 /* NDIMAX : format / dimension */
4632 /* NCF1MX : format by t of COURB1 */
4633 /* if ISENMSC= 1 : COURB1: The curve to be processed (NDIMAX,.) */
4634 /* NCOEFF : number of coeff of the curve */
4635 /* NCF2MX : format by t of COURB2 */
4636 /* NDIMEN : dimension of the curve and format of COURB2 */
4637 /* if ISENMSC=-1 : COURB2: The curve to be processed (.,NDIMEN) */
4639 /* OUTPUT ARGUMENTS : */
4640 /* --------------------- */
4641 /* if ISENMSC= 1 : COURB2: The resulting curve (.,NDIMEN) */
4642 /* if ISENMSC=-1 : COURB1: The resulting curve (NDIMAX,.) */
4644 /* COMMONS USED : */
4645 /* ------------------ */
4647 /* REFERENCES CALLED : */
4648 /* --------------------- */
4650 /* DESCRIPTION/NOTES/LIMITATIONS : */
4651 /* ----------------------------------- */
4652 /* allow to process the usual transfers as follows : */
4653 /* | ---- ISENMSC = 1 ---- | | ---- ISENMSC =-1 ----- | */
4654 /* TS (3,21) --> (21,3) AB ; AB (21,3) --> (3,21) TS */
4655 /* TS (3,21) --> (NU,3) AB ; AB (NU,3) --> (3,21) TS */
4656 /* (3,NU) --> (21,3) AB ; AB (21,3) --> (3,NU) */
4657 /* (3,NU) --> (NU,3) AB ; AB (NU,3) --> (3,NU) */
4659 /* ***********************************************************************
4663 /* Parameter adjustments */
4664 courb1_dim1 = *ndimax;
4665 courb1_offset = courb1_dim1 + 1;
4666 courb1 -= courb1_offset;
4667 courb2_dim1 = *ncf2mx;
4668 courb2_offset = courb2_dim1 + 1;
4669 courb2 -= courb2_offset;
4672 if (*ndimen > *ndimax || *ncoeff > *ncf1mx || *ncoeff > *ncf2mx) {
4676 if (*ndimen == 1 && *ncf1mx == *ncf2mx) {
4677 nboct = *ncf2mx << 3;
4678 if (*isenmsc == 1) {
4679 AdvApp2Var_SysBase::mcrfill_(&nboct,
4680 &courb1[courb1_offset],
4681 &courb2[courb2_offset]);
4683 if (*isenmsc == -1) {
4684 AdvApp2Var_SysBase::mcrfill_(&nboct,
4685 &courb2[courb2_offset],
4686 &courb1[courb1_offset]);
4693 if (*isenmsc == 1) {
4695 for (nd = 1; nd <= i__1; ++nd) {
4697 for (i__ = 1; i__ <= i__2; ++i__) {
4698 courb2[i__ + nd * courb2_dim1] = courb1[nd + i__ *
4704 } else if (*isenmsc == -1) {
4706 for (nd = 1; nd <= i__1; ++nd) {
4708 for (i__ = 1; i__ <= i__2; ++i__) {
4709 courb1[nd + i__ * courb1_dim1] = courb2[i__ + nd *
4721 /* ***********************************************************************
4729 AdvApp2Var_SysBase::maermsg_("MMFMCB5", iercod, 7L);
4734 //=======================================================================
4735 //function : AdvApp2Var_MathBase::mmfmtb1_
4737 //=======================================================================
4738 int AdvApp2Var_MathBase::mmfmtb1_(integer *maxsz1,
4750 /* System generated locals */
4751 integer table1_dim1, table1_offset, table2_dim1, table2_offset, i__1,
4754 /* Local variables */
4755 doublereal* work = 0;
4756 integer ilong, isize, ii, jj, ier = 0;
4757 intptr_t iofst = 0,iipt, jjpt;
4760 /************************************************************************
4765 /* Inversion of elements of a rectangular table (T1(i,j) */
4766 /* loaded in T2(j,i)) */
4770 /* ALL, MATH_ACCES :: TABLEAU&, INVERSION, &TABLEAU */
4772 /* INPUT ARGUMENTS : */
4773 /* ------------------ */
4774 /* MAXSZ1: Max Nb of elements by the 1st dimension of TABLE1. */
4775 /* TABLE1: Table of reals by two dimensions. */
4776 /* ISIZE1: Nb of useful elements of TABLE1 on the 1st dimension */
4777 /* JSIZE1: Nb of useful elements of TABLE1 on the 2nd dimension */
4778 /* MAXSZ2: Nb max of elements by the 1st dimension of TABLE2. */
4780 /* OUTPUT ARGUMENTS : */
4781 /* ------------------- */
4782 /* TABLE2: Table of reals by two dimensions, containing the transposition */
4783 /* of the rectangular table TABLE1. */
4784 /* ISIZE2: Nb of useful elements of TABLE2 on the 1st dimension */
4785 /* JSIZE2: Nb of useful elements of TABLE2 on the 2nd dimension */
4786 /* IERCOD: Erroe coder. */
4788 /* = 1, error in the dimension of tables */
4789 /* ether MAXSZ1 < ISIZE1 (table TABLE1 too small). */
4790 /* or MAXSZ2 < JSIZE1 (table TABLE2 too small). */
4792 /* COMMONS USED : */
4793 /* ---------------- */
4795 /* REFERENCES CALLED : */
4796 /* ---------------------- */
4798 /* DESCRIPTION/NOTES/LIMITATIONS : */
4799 /* ----------------------------------- */
4800 /* It is possible to use TABLE1 as input and output table i.e. */
4802 /* CALL MMFMTB1(MAXSZ1,TABLE1,ISIZE1,JSIZE1,MAXSZ2,TABLE1 */
4803 /* ,ISIZE2,JSIZE2,IERCOD) */
4806 /* **********************************************************************
4810 /* Parameter adjustments */
4811 table1_dim1 = *maxsz1;
4812 table1_offset = table1_dim1 + 1;
4813 table1 -= table1_offset;
4814 table2_dim1 = *maxsz2;
4815 table2_offset = table2_dim1 + 1;
4816 table2 -= table2_offset;
4817 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
4821 if (*isize1 > *maxsz1 || *jsize1 > *maxsz2) {
4826 isize = *maxsz2 * *isize1;
4827 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &isize, work, &iofst, &ier);
4832 /* DO NOT BE AFRAID OF CRUSHING. */
4835 for (ii = 1; ii <= i__1; ++ii) {
4836 iipt = (ii - 1) * *maxsz2 + iofst;
4838 for (jj = 1; jj <= i__2; ++jj) {
4839 jjpt = iipt + (jj - 1);
4840 work[jjpt] = table1[ii + jj * table1_dim1];
4846 AdvApp2Var_SysBase::mcrfill_(&ilong,
4848 &table2[table2_offset]);
4850 /* -------------- The number of elements of TABLE2 is returned ------------
4859 /* ------------------------------- THE END ------------------------------
4861 /* --> Invalid input. */
4865 /* --> Pb of allocation. */
4872 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &isize, work, &iofst, &ier);
4880 //=======================================================================
4881 //function : AdvApp2Var_MathBase::mmgaus1_
4883 //=======================================================================
4884 int AdvApp2Var_MathBase::mmgaus1_(integer *ndimf,
4901 /* System generated locals */
4904 /* Local variables */
4908 doublereal t, u[20], x;
4910 doublereal c1x, c2x;
4911 /* **********************************************************************
4917 /* Calculate the integral of function BFUNX passed in parameter */
4918 /* between limits XD and XF . */
4919 /* The function should be calculated for any value */
4920 /* of the variable in the given interval.. */
4921 /* The method GAUSS-LEGENDRE is used. */
4922 /* For explications refer to the book : */
4923 /* Complements de mathematiques a l'usage des Ingenieurs de */
4924 /* l'electrotechnique et des telecommunications. */
4925 /* Par Andre ANGOT - Collection technique et scientifique du CNET
4928 /* The degree of LEGENDRE polynoms used is passed in parameter.
4932 /* INTEGRATION,LEGENDRE,GAUSS */
4934 /* INPUT ARGUMENTS : */
4935 /* ------------------ */
4937 /* NDIMF : Dimension of the function */
4938 /* BFUNX : Function to integrate passed as argument */
4939 /* Should be declared as EXTERNAL in the call routine. */
4940 /* SUBROUTINE BFUNX(NDIMF,X,VAL,IER) */
4942 /* K : Parameter determining the degree of the LEGENDRE polynom that
4944 /* can take a value between 0 and 10. */
4945 /* The degree of the polynom is equal to 4 k, that is 4, 8,
4947 /* 12, 16, 20, 24, 28, 32, 36 and 40. */
4948 /* If K is not correct, the degree is set to 40 directly.
4950 /* XD : Lower limit of the interval of integration. */
4951 /* XF : Upper limit of the interval of integration. */
4952 /* SAUX1 : Auxiliary table */
4953 /* SAUX2 : Auxiliary table */
4955 /* OUTPUT ARGUMENTS : */
4956 /* ------------------- */
4958 /* SOMME : Value of the integral */
4959 /* NITER : Number of iterations to be carried out. */
4960 /* It is equal to the degree of the polynom. */
4962 /* IER : Error code : */
4963 /* < 0 ==> Attention - Warning */
4964 /* = 0 ==> Everything is OK */
4965 /* > 0 ==> Critical error - Apply special processing */
4966 /* ==> Error in the calculation of BFUNX (return code */
4967 /* of this routine */
4969 /* If error => SUM = 0 */
4971 /* COMMONS USED : */
4972 /* ----------------- */
4976 /* REFERENCES CALLED : */
4977 /* ---------------------- */
4980 /* @ BFUNX MVGAUS0 */
4982 /* DESCRIPTION/NOTES/LIMITATIONS : */
4983 /* --------------------------------- */
4985 /* See the explanations detailed in the listing */
4986 /* Use of the GAUSS method (orthogonal polynoms) */
4987 /* The symmetry of roots of these polynomes is used */
4988 /* Depending on K, the degree of the interpolated polynom grows.
4990 /* If you wish to calculate the integral with a given precision, */
4991 /* loop on k varying from 1 to 10 and test the difference of 2
4993 /* consecutive iterations. Stop the loop if this difference is less that */
4994 /* an epsilon value set to 10E-6 for example. */
4995 /* If S1 and S2 are 2 successive iterations, test following this example :
4998 /* AF=DABS(S1-S2) */
5000 /* If AS < 1 test if FS < eps otherwise test if AF/AS < eps
5002 /* -- ----- ----- */
5004 /************************************************************************
5007 /************************************************************************
5012 /* ****** General Initialization */
5014 /* Parameter adjustments */
5020 AdvApp2Var_SysBase::mvriraz_(ndimf,
5024 /* ****** Loading of coefficients U and H ** */
5025 /* -------------------------------------------- */
5027 mvgaus0_(k, u, h__, &ndeg, iercod);
5032 /* ****** C1X => Medium interval point [XD,XF] */
5033 /* ****** C2X => 1/2 amplitude interval [XD,XF] */
5035 c1x = (*xf + *xd) * .5;
5036 c2x = (*xf - *xd) * .5;
5038 /* ---------------------------------------- */
5039 /* ****** Integration for degree NDEG ** */
5040 /* ---------------------------------------- */
5043 for (j = 1; j <= i__1; ++j) {
5047 (*bfunx)(ndimf, &x, &saux1[1], iercod);
5053 (*bfunx)(ndimf, &x, &saux2[1], iercod);
5059 for (idimf = 1; idimf <= i__2; ++idimf) {
5060 somme[idimf] += h__[j - 1] * (saux1[idimf] + saux2[idimf]);
5067 for (idimf = 1; idimf <= i__1; ++idimf) {
5068 somme[idimf] *= c2x;
5071 /* ****** End of sub-program ** */
5077 //=======================================================================
5078 //function : mmherm0_
5080 //=======================================================================
5081 int mmherm0_(doublereal *debfin,
5084 integer c__576 = 576;
5088 /* System generated locals */
5092 /* Local variables */
5093 doublereal amat[36] /* was [6][6] */;
5096 integer iord1, iord2;
5097 doublereal miden[36] /* was [6][6] */;
5099 doublereal epspi, d1, d2;
5100 integer ii, jj, pp, ncf;
5102 integer iof[2], ier;
5103 doublereal mat[36] /* was [6][6] */;
5105 doublereal abid[72] /* was [12][6] */;
5106 /* ***********************************************************************
5111 /* INIT OF COEFFS. OF POLYNOMS OF HERMIT INTERPOLATION */
5115 /* MATH_ACCES :: HERMITE */
5117 /* INPUT ARGUMENTS */
5118 /* -------------------- */
5119 /* DEBFIN : PARAMETERS DEFINING THE CONSTRAINTS */
5120 /* DEBFIN(1) : FIRST PARAMETER */
5121 /* DEBFIN(2) : SECOND PARAMETER */
5123 /* ONE SHOULD HAVE: */
5124 /* ABS (DEBFIN(I)) < 100 */
5126 /* (ABS(DEBFIN(1)+ABS(DEBFIN(2))) > 1/100 */
5127 /* (for overflows) */
5129 /* ABS(DEBFIN(2)-DEBFIN(1)) / (ABS(DEBFIN(1)+ABS(DEBFIN(2))) > 1/100
5131 /* (for the conditioning) */
5134 /* OUTPUT ARGUMENTS : */
5135 /* --------------------- */
5137 /* IERCOD : Error code : 0 : O.K. */
5138 /* 1 : value of DEBFIN */
5139 /* are unreasonable */
5140 /* -1 : init was already done */
5141 /* (OK but no processing) */
5143 /* COMMONS USED : */
5144 /* ------------------ */
5146 /* REFERENCES CALLED : */
5147 /* ---------------------- */
5150 /* DESCRIPTION/NOTES/LIMITATIONS : */
5151 /* ----------------------------------- */
5153 /* This program initializes the coefficients of Hermit polynoms */
5154 /* that are read later by MMHERM1 */
5155 /* ***********************************************************************
5160 /* **********************************************************************
5165 /* Used to STORE coefficients of Hermit interpolation polynoms */
5171 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5172 /* ----------------------------------- */
5174 /* The coefficients of hermit polynoms are calculated by */
5175 /* the routine MMHERM0 and read by the routine MMHERM1 */
5177 /* **********************************************************************
5184 /* NBCOEF is the size of CMHERM (see below) */
5185 /* ***********************************************************************
5194 /* ***********************************************************************
5197 /* ***********************************************************************
5201 /* Parameter adjustments */
5205 d1 = advapp_abs(debfin[1]);
5206 if (d1 > (float)100.) {
5210 d2 = advapp_abs(debfin[2]);
5211 if (d2 > (float)100.) {
5216 if (d2 < (float).01) {
5220 d1 = (d__1 = debfin[2] - debfin[1], advapp_abs(d__1));
5221 if (d1 / d2 < (float).01) {
5226 /* ***********************************************************************
5228 /* Initialization */
5229 /* ***********************************************************************
5237 /* ***********************************************************************
5240 /* IS IT ALREADY INITIALIZED ? */
5242 d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
5245 if (debfin[1] != mmcmher_.tdebut) {
5248 if (debfin[2] != mmcmher_.tfinal) {
5251 if (d1 != mmcmher_.verifi) {
5259 /* ***********************************************************************
5262 /* ***********************************************************************
5268 /* Init. matrix identity : */
5271 AdvApp2Var_SysBase::mvriraz_(&ncmat,
5274 for (ii = 1; ii <= 6; ++ii) {
5275 miden[ii + ii * 6 - 7] = 1.;
5281 /* Init to 0 of table CMHERM */
5283 AdvApp2Var_SysBase::mvriraz_(&c__576, mmcmher_.cmherm);
5285 /* Calculation by solution of linear systems */
5287 for (iord1 = -1; iord1 <= 2; ++iord1) {
5288 for (iord2 = -1; iord2 <= 2; ++iord2) {
5295 iof[1] = iord[0] + 1;
5298 ncf = iord[0] + iord[1] + 2;
5300 /* Calculate matrix MAT to invert: */
5302 for (cot = 1; cot <= 2; ++cot) {
5305 if (iord[cot - 1] > -1) {
5308 for (jj = 1; jj <= i__1; ++jj) {
5314 i__1 = iord[cot - 1] + 1;
5315 for (pp = 1; pp <= i__1; ++pp) {
5317 ii = pp + iof[cot - 1];
5322 for (jj = 1; jj <= i__2; ++jj) {
5323 mat[ii + jj * 6 - 7] = (float)0.;
5328 for (jj = pp; jj <= i__2; ++jj) {
5330 /* everything is done in these 3 lines
5333 mat[ii + jj * 6 - 7] = cof[jj - 1] * prod;
5334 cof[jj - 1] *= jj - pp;
5335 prod *= debfin[cot];
5348 AdvApp2Var_MathBase::mmmrslwd_(&c__6, &ncf, &ncf, mat, miden, &epspi, abid, amat, &
5355 for (cot = 1; cot <= 2; ++cot) {
5356 i__1 = iord[cot - 1] + 1;
5357 for (pp = 1; pp <= i__1; ++pp) {
5359 for (ii = 1; ii <= i__2; ++ii) {
5360 mmcmher_.cmherm[ii + (pp + (cot + ((iord1 + (iord2 <<
5361 2)) << 1)) * 3) * 6 + 155] = amat[ii + (pp +
5362 iof[cot - 1]) * 6 - 7];
5375 /* ***********************************************************************
5378 /* The initialized flag is located: */
5380 mmcmher_.tdebut = debfin[1];
5381 mmcmher_.tfinal = debfin[2];
5383 d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
5384 mmcmher_.verifi = d1 * 16111959;
5387 /* ***********************************************************************
5392 /* ***********************************************************************
5403 /* ***********************************************************************
5408 AdvApp2Var_SysBase::maermsg_("MMHERM0", iercod, 7L);
5410 /* ***********************************************************************
5415 //=======================================================================
5416 //function : mmherm1_
5418 //=======================================================================
5419 int mmherm1_(doublereal *debfin,
5425 /* System generated locals */
5426 integer hermit_dim1, hermit_dim2, hermit_offset;
5428 /* Local variables */
5433 /* ***********************************************************************
5438 /* reading of coeffs. of HERMIT interpolation polynoms */
5442 /* MATH_ACCES :: HERMIT */
5444 /* INPUT ARGUMENTS : */
5445 /* -------------------- */
5446 /* DEBFIN : PARAMETES DEFINING THE CONSTRAINTS */
5447 /* DEBFIN(1) : FIRST PARAMETER */
5448 /* DEBFIN(2) : SECOND PARAMETER */
5450 /* Should be equal to the corresponding arguments during the */
5451 /* last call to MMHERM0 for the initialization of coeffs. */
5453 /* ORDRMX : indicates the dimensioning of HERMIT: */
5454 /* there is no choice : ORDRMX should be equal to the value */
5455 /* of PARAMETER IORDMX of INCLUDE MMCMHER, or 2 for the moment */
5457 /* IORDRE (2) : Orders of constraints in each corresponding parameter DEBFIN(I) */
5458 /* should be between -1 (no constraints) and ORDRMX. */
5461 /* OUTPUT ARGUMENTS : */
5462 /* --------------------- */
5464 /* HERMIT : HERMIT(1:IORDRE(1)+IORDRE(2)+2, j, cote) are the */
5465 /* coefficients in the canonic base of Hermit polynom */
5466 /* corresponding to orders IORDRE with parameters DEBFIN for */
5467 /* the constraint of order j on DEBFIN(cote). j is between 0 and IORDRE(cote). */
5470 /* IERCOD : Error code : */
5471 /* -1: O.K but necessary to reinitialize the coefficients */
5472 /* (info for optimization) */
5474 /* 1 : Error in MMHERM0 */
5475 /* 2 : arguments invalid */
5477 /* COMMONS USED : */
5478 /* ------------------ */
5480 /* REFERENCES CALLED : */
5481 /* ---------------------- */
5484 /* DESCRIPTION/NOTES/LIMITATIONS : */
5485 /* ----------------------------------- */
5487 /* This program reads coefficients of Hermit polynoms */
5488 /* that were earlier initialized by MMHERM0 */
5490 /* PMN : initialisation is no more done by the caller. */
5493 /* ***********************************************************************
5498 /* **********************************************************************
5503 /* Serves to STORE the coefficients of Hermit interpolation polynoms */
5509 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5510 /* ----------------------------------- */
5512 /* the coefficients of Hetmit polynoms are calculated by */
5513 /* routine MMHERM0 and read by routine MMHERM1 */
5516 /* **********************************************************************
5523 /* NBCOEF is the size of CMHERM (see lower) */
5527 /* ***********************************************************************
5534 /* ***********************************************************************
5536 /* Initializations */
5537 /* ***********************************************************************
5540 /* Parameter adjustments */
5542 hermit_dim1 = (*ordrmx << 1) + 2;
5543 hermit_dim2 = *ordrmx + 1;
5544 hermit_offset = hermit_dim1 * hermit_dim2 + 1;
5545 hermit -= hermit_offset;
5552 /* ***********************************************************************
5555 /* ***********************************************************************
5563 for (cot = 1; cot <= 2; ++cot) {
5564 if (iordre[cot] < -1) {
5567 if (iordre[cot] > *ordrmx) {
5574 /* IS-IT CORRECTLY INITIALIZED ? */
5576 d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
5579 /* OTHERWISE IT IS INITIALIZED */
5581 if (debfin[1] != mmcmher_.tdebut || debfin[2] != mmcmher_.tfinal || d1
5582 != mmcmher_.verifi) {
5584 mmherm0_(&debfin[1], iercod);
5591 /* ***********************************************************************
5594 /* ***********************************************************************
5599 AdvApp2Var_SysBase::msrfill_(&nbval, &mmcmher_.cmherm[((((iordre[1] + (iordre[2] << 2)) << 1)
5600 + 1) * 3 + 1) * 6 + 156], &hermit[hermit_offset]);
5602 /* ***********************************************************************
5607 /* ***********************************************************************
5618 /* ***********************************************************************
5623 AdvApp2Var_SysBase::maermsg_("MMHERM1", iercod, 7L);
5625 /* ***********************************************************************
5630 //=======================================================================
5631 //function : AdvApp2Var_MathBase::mmhjcan_
5633 //=======================================================================
5634 int AdvApp2Var_MathBase::mmhjcan_(integer *ndimen,
5647 /* System generated locals */
5648 integer tcbold_dim1, tcbold_dim2, tcbold_offset, tcbnew_dim1, tcbnew_dim2,
5649 tcbnew_offset, i__1, i__2, i__3, i__4, i__5;
5652 /* Local variables */
5655 doublereal taux1[21];
5656 integer d__, e, i__, k;
5659 doublereal tjacap[21];
5661 doublereal hermit[36]/* was [6][3][2] */, ctenor, bornes[2];
5665 /* ***********************************************************************
5670 /* CONVERSION OF TABLE TCBOLD OF POLYNOMIAL CURVE COEFFICIENTS */
5671 /* EXPRESSED IN HERMIT JACOBI BASE, INTO A */
5672 /* TABLE OF COEFFICIENTS TCBNEW OF COURVES EXPRESSED IN THE CANONIC BASE */
5676 /* CANNONIC, HERMIT, JACCOBI */
5678 /* INPUT ARGUMENTS : */
5679 /* -------------------- */
5680 /* ORDHER : ORDER OF HERMIT POLYNOMS OR ORDER OF CONTINUITY */
5681 /* NCOEFS : NUMBER OF COEFFICIENTS OF A POLYNOMIAL CURVE */
5682 /* FOR ONE OF ITS NDIM COMPONENTS;(DEGREE+1 OF THE CURVE)
5684 /* NDIM : DIMENSION OF THE CURVE */
5685 /* CBHEJA : TABLE OF COEFFICIENTS OF THE CURVE IN THE BASE */
5687 /* (H(0,-1),..,H(ORDHER,-1),H(0,1),..,H(ORDHER,1), */
5688 /* JA(ORDHER+1,2*ORDHER+2),....,JA(ORDHER+1,NCOEFS-1) */
5690 /* OUTPUT ARGUMENTS : */
5691 /* --------------------- */
5692 /* CBRCAN : TABLE OF COEFFICIENTS OF THE CURVE IN THE CANONIC BASE */
5695 /* COMMONS USED : */
5696 /* ------------------ */
5699 /* REFERENCES CALLED : */
5700 /* --------------------- */
5703 /* ***********************************************************************
5707 /* ***********************************************************************
5712 /* Providesinteger constants from 0 to 1000 */
5718 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5719 /* ----------------------------------- */
5721 /* ***********************************************************************
5725 /* ***********************************************************************
5731 /* ***********************************************************************
5733 /* INITIALIZATION */
5734 /* ***********************************************************************
5737 /* Parameter adjustments */
5739 tcbnew_dim1 = *ndimen;
5740 tcbnew_dim2 = *ncflim;
5741 tcbnew_offset = tcbnew_dim1 * (tcbnew_dim2 + 1) + 1;
5742 tcbnew -= tcbnew_offset;
5743 tcbold_dim1 = *ndimen;
5744 tcbold_dim2 = *ncflim;
5745 tcbold_offset = tcbold_dim1 * (tcbold_dim2 + 1) + 1;
5746 tcbold -= tcbold_offset;
5749 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
5751 AdvApp2Var_SysBase::mgenmsg_("MMHJCAN", 7L);
5758 /* ***********************************************************************
5761 /* ***********************************************************************
5771 /* CALCULATION OF HERMIT POLYNOMS IN THE CANONIC BASE ON (-1,1) */
5774 iordre[0] = *orcont;
5775 iordre[1] = *orcont;
5776 mmherm1_(bornes, &c__2, iordre, hermit, &ier);
5786 for (e = 1; e <= i__1; ++e) {
5788 ctenor = (tdecop[e] - tdecop[e - 1]) / 2;
5796 for (d__ = 1; d__ <= i__2; ++d__) {
5798 /* CONVERSION OF THE COEFFICIENTS OF THE PART OF THE CURVE EXPRESSED */
5799 /* IN HERMIT BASE, INTO THE CANONIC BASE */
5801 AdvApp2Var_SysBase::mvriraz_(&ncoeff, taux1);
5804 for (k = 1; k <= i__3; ++k) {
5806 for (i__ = 1; i__ <= i__4; ++i__) {
5808 mfact = AdvApp2Var_MathBase::pow__di(&ctenor, &i__5);
5809 taux1[k - 1] += (tcbold[d__ + (i__ + e * tcbold_dim2) *
5810 tcbold_dim1] * hermit[k + (i__ + 2) * 6 - 19] +
5811 tcbold[d__ + (i__ + aux1 + e * tcbold_dim2) *
5812 tcbold_dim1] * hermit[k + (i__ + 5) * 6 - 19]) *
5819 for (i__ = aux2 + 1; i__ <= i__3; ++i__) {
5820 taux1[i__ - 1] = tcbold[d__ + (i__ + e * tcbold_dim2) *
5824 /* CONVERSION OF THE COEFFICIENTS OF THE PART OF THE CURVE EXPRESSED */
5825 /* IN CANONIC-JACOBI BASE, INTO THE CANONIC BASE */
5829 AdvApp2Var_MathBase::mmapcmp_(&minombr_.nbr[1], &c__21, &ncoeff, taux1, tjacap);
5830 AdvApp2Var_MathBase::mmjacan_(orcont, &ndeg, tjacap, taux1);
5832 /* RECOPY THE COEFS RESULTING FROM THE CONVERSION IN THE TABLE */
5836 for (i__ = 1; i__ <= i__3; ++i__) {
5837 tcbnew[d__ + (i__ + e * tcbnew_dim2) * tcbnew_dim1] = taux1[
5846 /* ***********************************************************************
5848 /* PROCESSING OF ERRORS */
5849 /* ***********************************************************************
5859 /* ***********************************************************************
5861 /* RETURN CALLING PROGRAM */
5862 /* ***********************************************************************
5867 AdvApp2Var_SysBase::maermsg_("MMHJCAN", iercod, 7L);
5869 AdvApp2Var_SysBase::mgsomsg_("MMHJCAN", 7L);
5874 //=======================================================================
5875 //function : AdvApp2Var_MathBase::mminltt_
5877 //=======================================================================
5878 int AdvApp2Var_MathBase::mminltt_(integer *ncolmx,
5884 doublereal *,//epseg,
5887 /* System generated locals */
5888 integer tabtri_dim1, tabtri_offset, i__1, i__2;
5890 /* Local variables */
5892 integer icol, ilgn, nlgn, noct, inser;
5893 doublereal epsega = 0.;
5896 /* ***********************************************************************
5901 /* . Insert a line in a table parsed without redundance */
5905 /* TOUS,MATH_ACCES :: TABLEAU&,INSERTION,&TABLEAU */
5907 /* INPUT ARGUMENTS : */
5908 /* -------------------- */
5909 /* . NCOLMX : Number of columns in the table */
5910 /* . NLGNMX : Number of lines in the table */
5911 /* . TABTRI : Table parsed by lines without redundances */
5912 /* . NBRCOL : Number of columns used */
5913 /* . NBRLGN : Number of lines used */
5914 /* . AJOUTE : Line to be added */
5915 /* . EPSEGA : Epsilon to test the redundance */
5917 /* OUTPUT ARGUMENTS : */
5918 /* --------------------- */
5919 /* . TABTRI : Table parsed by lines without redundances */
5920 /* . NBRLGN : Number of lines used */
5921 /* . IERCOD : 0 -> No problem */
5922 /* 1 -> The table is full */
5924 /* COMMONS USED : */
5925 /* ------------------ */
5927 /* REFERENCES CALLED : */
5928 /* --------------------- */
5930 /* DESCRIPTION/NOTES/LIMITATIONS : */
5931 /* ----------------------------------- */
5932 /* . The line is inserted only if there is no line with all
5934 /* elements equl to those which are planned to be insered, to epsilon. */
5936 /* . Level of de debug = 3 */
5940 /* DECLARATIONS , CONTROL OF INPUT ARGUMENTS , INITIALIZATION */
5941 /* ***********************************************************************
5944 /* --- Parameters */
5950 /* --- Local variables */
5955 /* Parameter adjustments */
5956 tabtri_dim1 = *ncolmx;
5957 tabtri_offset = tabtri_dim1 + 1;
5958 tabtri -= tabtri_offset;
5962 ibb = AdvApp2Var_SysBase::mnfndeb_();
5965 AdvApp2Var_SysBase::mgenmsg_("MMINLTT", 7L);
5968 /* --- Control arguments */
5970 if (*nbrlgn >= *nlgnmx) {
5974 /* -------------------- */
5975 /* *** INITIALIZATION */
5976 /* -------------------- */
5980 /* ---------------------------- */
5981 /* *** SEARCH OF REDUNDANCE */
5982 /* ---------------------------- */
5985 for (ilgn = 1; ilgn <= i__1; ++ilgn) {
5986 if (tabtri[ilgn * tabtri_dim1 + 1] >= ajoute[1] - epsega) {
5987 if (tabtri[ilgn * tabtri_dim1 + 1] <= ajoute[1] + epsega) {
5989 for (icol = 1; icol <= i__2; ++icol) {
5990 if (tabtri[icol + ilgn * tabtri_dim1] < ajoute[icol] -
5991 epsega || tabtri[icol + ilgn * tabtri_dim1] >
5992 ajoute[icol] + epsega) {
6006 /* ----------------------------------- */
6007 /* *** SEARCH OF THE INSERTION POINT */
6008 /* ----------------------------------- */
6013 for (ilgn = 1; ilgn <= i__1; ++ilgn) {
6015 for (icol = 1; icol <= i__2; ++icol) {
6016 if (tabtri[icol + ilgn * tabtri_dim1] < ajoute[icol]) {
6019 if (tabtri[icol + ilgn * tabtri_dim1] > ajoute[icol]) {
6030 /* -------------- */
6032 /* -------------- */
6039 /* --- Shift lower */
6041 nlgn = *nbrlgn - inser;
6043 noct = (*ncolmx << 3) * nlgn;
6044 AdvApp2Var_SysBase::mcrfill_(&noct,
6045 &tabtri[inser * tabtri_dim1 + 1],
6046 &tabtri[(inser + 1)* tabtri_dim1 + 1]);
6051 noct = *nbrcol << 3;
6052 AdvApp2Var_SysBase::mcrfill_(&noct,
6054 &tabtri[inser * tabtri_dim1 + 1]);
6058 /* ******************************************************************** */
6059 /* OUTPUT ERROR , RETURN CALLING PROGRAM , MESSAGES */
6060 /* ******************************************************************** */
6062 /* --- The table is already full */
6071 AdvApp2Var_SysBase::maermsg_("MMINLTT", iercod, 7L);
6074 AdvApp2Var_SysBase::mgsomsg_("MMINLTT", 7L);
6079 //=======================================================================
6080 //function : AdvApp2Var_MathBase::mmjacan_
6082 //=======================================================================
6083 int AdvApp2Var_MathBase::mmjacan_(const integer *ideriv,
6088 /* System generated locals */
6089 integer poljac_dim1, i__1, i__2;
6091 /* Local variables */
6092 integer iptt, i__, j, ibb;
6095 /* ***********************************************************************
6100 /* Routine of transfer of Jacobi normalized to canonic [-1,1], */
6101 /* the tables are ranked by even, then by uneven degree. */
6105 /* LEGENDRE,JACOBI,PASSAGE. */
6107 /* INPUT ARGUMENTS : */
6108 /* ------------------ */
6109 /* IDERIV : Order of Jacobi between -1 and 2. */
6110 /* NDEG : The true degree of the polynom. */
6111 /* POLJAC : The polynom in the Jacobi base. */
6113 /* OUTPUT ARGUMENTS : */
6114 /* ------------------- */
6115 /* POLCAN : The curve expressed in the canonic base [-1,1]. */
6117 /* COMMONS USED : */
6118 /* ---------------- */
6120 /* REFERENCES CALLED : */
6121 /* ----------------------- */
6123 /* DESCRIPTION/NOTES/LIMITATIONS : */
6124 /* ----------------------------------- */
6127 /* ***********************************************************************
6130 /* Name of the routine */
6132 /* Matrices of conversion */
6135 /* ***********************************************************************
6140 /* MATRIX OF TRANSFORMATION OF LEGENDRE BASE */
6146 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
6147 /* ----------------------------------- */
6150 /* ***********************************************************************
6155 /* Legendre common / Restricted Casteljau. */
6157 /* 0:1 0 Concerns the even terms, 1 the uneven terms. */
6158 /* CANPLG : Matrix of passage to canonic from Jacobi with calculated parities */
6159 /* PLGCAN : Matrix of passage from Jacobi to canonic with calculated parities */
6162 /* ***********************************************************************
6165 /* Parameter adjustments */
6166 poljac_dim1 = *ndeg / 2 + 1;
6169 ibb = AdvApp2Var_SysBase::mnfndeb_();
6171 AdvApp2Var_SysBase::mgenmsg_("MMJACAN", 7L);
6174 /* ----------------- Expression of terms of even degree ----------------
6178 for (i__ = 0; i__ <= i__1; ++i__) {
6180 iptt = i__ * 31 - (i__ + 1) * i__ / 2 + 1;
6182 for (j = i__; j <= i__2; ++j) {
6183 bid += mmjcobi_.plgcan[iptt + j + *ideriv * 992 + 991] * poljac[
6187 polcan[i__ * 2] = bid;
6191 /* --------------- Expression of terms of uneven degree ----------------
6198 i__1 = (*ndeg - 1) / 2;
6199 for (i__ = 0; i__ <= i__1; ++i__) {
6201 iptt = i__ * 31 - (i__ + 1) * i__ / 2 + 1;
6202 i__2 = (*ndeg - 1) / 2;
6203 for (j = i__; j <= i__2; ++j) {
6204 bid += mmjcobi_.plgcan[iptt + j + ((*ideriv << 1) + 1) * 496 +
6205 991] * poljac[j + poljac_dim1];
6208 polcan[(i__ << 1) + 1] = bid;
6212 /* -------------------------------- The end -----------------------------
6217 AdvApp2Var_SysBase::mgsomsg_("MMJACAN", 7L);
6222 //=======================================================================
6223 //function : AdvApp2Var_MathBase::mmjaccv_
6225 //=======================================================================
6226 int AdvApp2Var_MathBase::mmjaccv_(const integer *ncoef,
6227 const integer *ndim,
6228 const integer *ider,
6229 const doublereal *crvlgd,
6234 /* Initialized data */
6236 static char nomprg[8+1] = "MMJACCV ";
6238 /* System generated locals */
6239 integer crvlgd_dim1, crvlgd_offset, crvcan_dim1, crvcan_offset,
6240 polaux_dim1, i__1, i__2;
6242 /* Local variables */
6243 integer ndeg, i__, nd, ii, ibb;
6245 /* ***********************************************************************
6250 /* Passage from the normalized Jacobi base to the canonic base. */
6254 /* SMOOTHING, BASE, LEGENDRE */
6257 /* INPUT ARGUMENTS : */
6258 /* ------------------ */
6259 /* NDIM: Space Dimension. */
6260 /* NCOEF: Degree +1 of the polynom. */
6261 /* IDER: Order of Jacobi polynoms. */
6262 /* CRVLGD : Curve in the base of Jacobi. */
6264 /* OUTPUT ARGUMENTS : */
6265 /* ------------------- */
6266 /* POLAUX : Auxilliary space. */
6267 /* CRVCAN : The curve in the canonic base [-1,1] */
6269 /* COMMONS USED : */
6270 /* ---------------- */
6272 /* REFERENCES CALLED : */
6273 /* ----------------------- */
6275 /* DESCRIPTION/NOTES/LIMITATIONS : */
6276 /* ----------------------------------- */
6279 /* *********************************************************************
6282 /* Name of the routine */
6283 /* Parameter adjustments */
6284 polaux_dim1 = (*ncoef - 1) / 2 + 1;
6285 crvcan_dim1 = *ncoef - 1 + 1;
6286 crvcan_offset = crvcan_dim1;
6287 crvcan -= crvcan_offset;
6288 crvlgd_dim1 = *ncoef - 1 + 1;
6289 crvlgd_offset = crvlgd_dim1;
6290 crvlgd -= crvlgd_offset;
6294 ibb = AdvApp2Var_SysBase::mnfndeb_();
6296 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
6302 for (nd = 1; nd <= i__1; ++nd) {
6303 /* Loading of the auxilliary table. */
6306 for (i__ = 0; i__ <= i__2; ++i__) {
6307 polaux[i__] = crvlgd[ii + nd * crvlgd_dim1];
6314 i__2 = (ndeg - 1) / 2;
6315 for (i__ = 0; i__ <= i__2; ++i__) {
6316 polaux[i__ + polaux_dim1] = crvlgd[ii + nd * crvlgd_dim1];
6321 /* Call the routine of base change. */
6322 AdvApp2Var_MathBase::mmjacan_(ider, &ndeg, polaux, &crvcan[nd * crvcan_dim1]);
6331 //=======================================================================
6332 //function : mmloncv_
6334 //=======================================================================
6335 int mmloncv_(integer *ndimax,
6345 /* Initialized data */
6349 /* System generated locals */
6350 integer courbe_dim1, courbe_offset, i__1, i__2;
6352 /* Local variables */
6355 doublereal c1, c2, d1, d2,
6356 wgaus[20] = {0.}, uroot[20] = {0.}, x1, x2, dd;
6359 doublereal der1, der2;
6364 /* **********************************************************************
6367 /* FUNCTION : Length of an arc of curve on a given interval */
6368 /* ---------- for a function the mathematic representation */
6369 /* which of is a multidimensional polynom. */
6370 /* The polynom is a set of polynoms the coefficients which of are ranked */
6371 /* in a table with 2 indices, each line relative to 1 polynom. */
6372 /* The polynom is defined by its coefficients ordered by increasing
6373 * power of the variable. */
6374 /* All polynoms have the same number of coefficients (and the same degree). */
6376 /* KEYWORDS : LENGTH, CURVE */
6379 /* INPUT ARGUMENTS : */
6380 /* -------------------- */
6382 /* NDIMAX : Max number of lines of tables (max number of polynoms). */
6383 /* NDIMEN : Dimension of the polynom (Nomber of polynoms). */
6384 /* NCOEFF : Number of coefficients of the polynom (no limitation) */
6385 /* This is degree + 1 */
6386 /* COURBE : Coefficients of the polynom ordered by increasing power */
6387 /* Dimension to (NDIMAX,NCOEFF). */
6388 /* TDEBUT : Lower limit of integration for length calculation. */
6389 /* TFINAL : Upper limit of integration for length calculation. */
6391 /* OUTPUT ARGUMENTS : */
6392 /* --------------------- */
6393 /* XLONGC : Length of arc of curve */
6395 /* IERCOD : Error code : */
6396 /* = 0 ==> All is OK */
6397 /* = 1 ==> NDIMEN or NCOEFF negative or null */
6398 /* = 2 ==> Pb loading Legendre roots and Gauss weight */
6401 /* If error => XLONGC = 0 */
6403 /* COMMONS USED : */
6404 /* ------------------ */
6408 /* REFERENCES CALLED : */
6409 /* ---------------------- */
6411 /* MAERMSG R*8 DSQRT I*4 MIN */
6414 /* DESCRIPTION/NOTES/LIMITATIONS : */
6415 /* ----------------------------------- */
6417 /* See VGAUSS to understand well the technique. */
6418 /* Actually SQRT (dpi^2) is integrated for i=1,nbdime */
6419 /* Calculation of the derivative is included in the code to avoid an additional */
6420 /* call of the routine. */
6422 /* The integrated function is strictly increasing, it */
6423 /* is not necessary to use a high degree for the GAUSS method GAUSS. */
6425 /* The degree of LEGENDRE polynom results from the degree of the */
6426 /* polynom to be integrated. It can vary from 4 to 40 (with step of 4). */
6428 /* The precision (relative) of integration is of order 1.D-8. */
6430 /* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
6432 /* Attention : the precision of the result is not controlled. */
6433 /* If you wish to control it, use MMCGLC1, taking into account that */
6434 /* the performance (in time) will be worse. */
6436 /* >=====================================================================
6439 /* ATTENTION : SAVE KGAR WGAUS and UROOT EVENTUALLY */
6441 /* INTEGER I1,I20 */
6442 /* PARAMETER (I1=1,I20=20) */
6444 /* Parameter adjustments */
6445 courbe_dim1 = *ndimax;
6446 courbe_offset = courbe_dim1 + 1;
6447 courbe -= courbe_offset;
6451 /* ****** General initialization ** */
6456 /* ****** Initialization of UROOT, WGAUS, NGAUS and KGAR ** */
6458 /* CALL MXVINIT(IERXV,'INTEGER',I1,KGAR,'INTEGER',I1,NGAUS */
6459 /* 1 ,'DOUBLE PRECISION',I20,UROOT,'DOUBLE PRECISION',I20,WGAUS) */
6460 /* IF (IERXV.GT.0) KGAR=0 */
6462 /* ****** Test the equity of limits ** */
6464 if (*tdebut == *tfinal) {
6469 /* ****** Test the dimension and the number of coefficients ** */
6471 if (*ndimen <= 0 || *ncoeff <= 0) {
6476 /* ****** Calculate the optimal degree ** */
6478 kk = *ncoeff / 4 + 1;
6479 kk = advapp_min(kk,10);
6481 /* ****** Return the coefficients for the integral (DEGRE=4*KK) */
6482 /* if KK <> KGAR. */
6485 mvgaus0_(&kk, uroot, wgaus, &ngaus, iercod);
6494 /* C1 => Point medium interval */
6495 /* C2 => 1/2 amplitude interval */
6497 c1 = (*tfinal + *tdebut) * .5;
6498 c2 = (*tfinal - *tdebut) * .5;
6500 /* ----------------------------------------------------------- */
6501 /* ****** Integration - Loop on GAUSS intervals ** */
6502 /* ----------------------------------------------------------- */
6507 for (jj = 1; jj <= i__1; ++jj) {
6509 /* ****** Integration taking the symmetry into account ** */
6511 tran = c2 * uroot[jj - 1];
6515 /* ****** Derivation on the dimension of the space ** */
6520 for (kk = 1; kk <= i__2; ++kk) {
6521 d1 = (*ncoeff - 1) * courbe[kk + *ncoeff * courbe_dim1];
6523 for (ii = *ncoeff - 1; ii >= 2; --ii) {
6524 dd = (ii - 1) * courbe[kk + ii * courbe_dim1];
6534 /* ****** Integration ** */
6536 som += wgaus[jj - 1] * c2 * (sqrt(der1) + sqrt(der2));
6538 /* ****** End of loop on GAUSS intervals ** */
6543 /* ****** Work ended ** */
6547 /* ****** It is forced IERCOD = 0 ** */
6551 /* ****** Final processing ** */
6555 /* ****** Save UROOT, WGAUS, NGAUS and KGAR ** */
6557 /* CALL MXVSAVE(IERXV,'INTEGER',I1,KGAR,'INTEGER',I1,NGAUS */
6558 /* 1 ,'DOUBLE PRECISION',I20,UROOT,'DOUBLE PRECISION',I20,WGAUS) */
6559 /* IF (IERXV.GT.0) KGAR=0 */
6561 /* ****** End of sub-program ** */
6564 AdvApp2Var_SysBase::maermsg_("MMLONCV", iercod, 7L);
6569 //=======================================================================
6570 //function : AdvApp2Var_MathBase::mmpobas_
6572 //=======================================================================
6573 int AdvApp2Var_MathBase::mmpobas_(doublereal *tparam,
6585 /* Initialized data */
6587 doublereal moin11[2] = { -1.,1. };
6589 /* System generated locals */
6590 integer valbas_dim1, i__1;
6592 /* Local variables */
6593 doublereal vjac[80], herm[24];
6596 integer nwcof, iunit;
6597 doublereal wpoly[7];
6598 integer ii, jj, iorjac;
6599 doublereal hermit[36] /* was [6][3][2] */;
6600 integer kk1, kk2, kk3;
6604 /* ***********************************************************************
6609 /* Position on the polynoms of base hermit-Jacobi */
6610 /* and their succesive derivatives */
6614 /* PUBLIC, POSITION, HERMIT, JACOBI */
6616 /* INPUT ARGUMENTS : */
6617 /* -------------------- */
6618 /* TPARAM : Parameter for which the position is found. */
6619 /* IORDRE : Orderof hermit-Jacobi (-1,0,1, ou 2) */
6620 /* NCOEFF : Number of coefficients of polynoms (Nb of value to calculate) */
6621 /* NDERIV : Number of derivative to calculate (0<= N <=3) */
6622 /* 0 -> Position simple on base functions */
6623 /* N -> Position on base functions and derivative */
6624 /* of order 1 to N */
6626 /* OUTPUT ARGUMENTS : */
6627 /* --------------------- */
6628 /* VALBAS (NCOEFF, 0:NDERIV) : calculated value */
6630 /* d vj(t) = VALBAS(J, I) */
6634 /* IERCOD : Error code */
6636 /* 1 : Incoherence of input arguments */
6638 /* COMMONS USED : */
6639 /* -------------- */
6642 /* REFERENCES CALLED : */
6643 /* ------------------- */
6646 /* DESCRIPTION/NOTES/LIMITATIONS : */
6647 /* ----------------------------------- */
6650 /* ***********************************************************************
6653 /* ***********************************************************************
6658 /* Parameter adjustments */
6659 valbas_dim1 = *ncoeff;
6664 /* ***********************************************************************
6666 /* INITIALIZATIONS */
6667 /* ***********************************************************************
6672 /* ***********************************************************************
6675 /* ***********************************************************************
6690 iorjac = (*iordre + 1) << 1;
6692 /* (1) Generic Calculations .... */
6694 /* (1.a) Calculation of hermit polynoms */
6697 mmherm1_(moin11, &c__2, iord, hermit, &ier);
6703 /* (1.b) Evaluation of hermit polynoms */
6706 iunit = *nderiv + 1;
6707 khe = (*iordre + 1) * iunit;
6712 for (ii = 0; ii <= i__1; ++ii) {
6713 mmdrvcb_(nderiv, &c__1, &iorjac, &hermit[(ii + 3) * 6 - 18],
6714 tparam, &herm[jj - 1], &ier);
6719 mmdrvcb_(nderiv, &c__1, &iorjac, &hermit[(ii + 6) * 6 - 18],
6720 tparam, &herm[jj + khe - 1], &ier);
6730 for (ii = 0; ii <= i__1; ++ii) {
6731 AdvApp2Var_MathBase::mmpocrb_(&c__1, &iorjac, &hermit[(ii + 3) * 6 - 18], &c__1,
6732 tparam, &herm[jj - 1]);
6734 AdvApp2Var_MathBase::mmpocrb_(&c__1, &iorjac, &hermit[(ii + 6) * 6 - 18], &c__1,
6735 tparam, &herm[jj + khe - 1]);
6740 /* (1.c) Evaluation of Jacobi polynoms */
6742 ii = *ncoeff - iorjac;
6744 mmpojac_(tparam, &iorjac, &ii, nderiv, vjac, &ier);
6749 /* (1.d) Evaluation of W(t) */
6753 nwcof = advapp_max(i__1,1);
6754 AdvApp2Var_SysBase::mvriraz_(&nwcof,
6761 } else if (*iordre == 1) {
6764 } else if (*iordre == 0) {
6768 mmdrvcb_(nderiv, &c__1, &nwcof, wpoly, tparam, wval, &ier);
6773 kk1 = *ncoeff - iorjac;
6777 /* (2) Evaluation of order 0 */
6781 for (ii = 1; ii <= i__1; ++ii) {
6782 valbas[ii] = herm[jj - 1];
6787 for (ii = 1; ii <= i__1; ++ii) {
6788 valbas[ii + iorjac] = wval[0] * vjac[ii - 1];
6791 /* (3) Evaluation of order 1 */
6796 for (ii = 1; ii <= i__1; ++ii) {
6797 valbas[ii + valbas_dim1] = herm[jj - 1];
6803 for (ii = 1; ii <= i__1; ++ii) {
6804 valbas[ii + iorjac + valbas_dim1] = wval[0] * vjac[ii + kk1 - 1]
6805 + wval[1] * vjac[ii - 1];
6809 /* (4) Evaluation of order 2 */
6814 for (ii = 1; ii <= i__1; ++ii) {
6815 valbas[ii + (valbas_dim1 << 1)] = herm[jj - 1];
6820 for (ii = 1; ii <= i__1; ++ii) {
6821 valbas[ii + iorjac + (valbas_dim1 << 1)] = wval[0] * vjac[ii +
6822 kk2 - 1] + wval[1] * 2 * vjac[ii + kk1 - 1] + wval[2] *
6827 /* (5) Evaluation of order 3 */
6832 for (ii = 1; ii <= i__1; ++ii) {
6833 valbas[ii + valbas_dim1 * 3] = herm[jj - 1];
6838 for (ii = 1; ii <= i__1; ++ii) {
6839 valbas[ii + iorjac + valbas_dim1 * 3] = wval[0] * vjac[ii + kk3 -
6840 1] + wval[1] * 3 * vjac[ii + kk2 - 1] + wval[2] * 3 *
6841 vjac[ii + kk1 - 1] + wval[3] * vjac[ii - 1];
6847 /* ***********************************************************************
6849 /* ERROR PROCESSING */
6850 /* ***********************************************************************
6860 /* ***********************************************************************
6862 /* RETURN CALLING PROGRAM */
6863 /* ***********************************************************************
6869 AdvApp2Var_SysBase::maermsg_("MMPOBAS", iercod, 7L);
6874 //=======================================================================
6875 //function : AdvApp2Var_MathBase::mmpocrb_
6877 //=======================================================================
6878 int AdvApp2Var_MathBase::mmpocrb_(integer *ndimax,
6886 /* System generated locals */
6887 integer courbe_dim1, courbe_offset, i__1, i__2;
6889 /* Local variables */
6891 integer isize, nd, kcf, ncf;
6894 /* ***********************************************************************
6899 /* CALCULATE THE COORDINATES OF A POINT OF A CURVE OF GIVEN PARAMETER */
6900 /* TPARAM ( IN 2D, 3D OR MORE) */
6904 /* TOUS , MATH_ACCES :: COURBE&,PARAMETRE& , POSITIONNEMENT , &POINT
6907 /* INPUT ARGUMENTS : */
6908 /* ------------------ */
6909 /* NDIMAX : format / dimension of the curve */
6910 /* NCOEFF : Nb of coefficients of the curve */
6911 /* COURBE : Matrix of coefficients of the curve */
6912 /* NDIM : Dimension useful of the workspace */
6913 /* TPARAM : Value of the parameter where the point is calculated */
6915 /* OUTPUT ARGUMENTS : */
6916 /* ------------------- */
6917 /* PNTCRB : Coordinates of the calculated point */
6919 /* COMMONS USED : */
6920 /* ---------------- */
6924 /* REFERENCES CALLED : */
6925 /* ---------------------- */
6927 /* MIRAZ MVPSCR2 MVPSCR3 */
6929 /* DESCRIPTION/NOTES/LIMITATIONS : */
6930 /* ----------------------------------- */
6933 /* ***********************************************************************
6937 /* ***********************************************************************
6940 /* Parameter adjustments */
6941 courbe_dim1 = *ndimax;
6942 courbe_offset = courbe_dim1 + 1;
6943 courbe -= courbe_offset;
6948 AdvApp2Var_SysBase::miraz_(&isize,
6955 /* optimal processing 3d */
6957 if (*ndim == 3 && *ndimax == 3) {
6958 mvpscr3_(ncoeff, &courbe[courbe_offset], tparam, &pntcrb[1]);
6960 /* optimal processing 2d */
6962 } else if (*ndim == 2 && *ndimax == 2) {
6963 mvpscr2_(ncoeff, &courbe[courbe_offset], tparam, &pntcrb[1]);
6965 /* Any dimension - scheme of HORNER */
6967 } else if (*tparam == 0.) {
6969 for (nd = 1; nd <= i__1; ++nd) {
6970 pntcrb[nd] = courbe[nd + courbe_dim1];
6973 } else if (*tparam == 1.) {
6975 for (ncf = 1; ncf <= i__1; ++ncf) {
6977 for (nd = 1; nd <= i__2; ++nd) {
6978 pntcrb[nd] += courbe[nd + ncf * courbe_dim1];
6984 ncof2 = *ncoeff + 2;
6986 for (nd = 1; nd <= i__1; ++nd) {
6988 for (ncf = 2; ncf <= i__2; ++ncf) {
6990 pntcrb[nd] = (pntcrb[nd] + courbe[nd + kcf * courbe_dim1]) * *
6994 pntcrb[nd] += courbe[nd + courbe_dim1];
7003 //=======================================================================
7004 //function : AdvApp2Var_MathBase::mmmpocur_
7006 //=======================================================================
7007 int AdvApp2Var_MathBase::mmmpocur_(integer *ncofmx,
7015 /* System generated locals */
7016 integer courbe_dim1, courbe_offset, i__1;
7018 /* Local variables */
7023 /* ***********************************************************************
7028 /* Position of a point on curve (ncofmx,ndim). */
7032 /* TOUS , AB_SPECIFI :: COURBE&,POLYNOME&,POSITIONNEMENT,&POINT */
7034 /* INPUT ARGUMENTS : */
7035 /* ------------------ */
7036 /* NCOFMX: Format / degree of the CURVE. */
7037 /* NDIM : Dimension of the space. */
7038 /* NDEG : Degree of the polynom. */
7039 /* COURBE: Coefficients of the curve. */
7040 /* TPARAM: Parameter on the curve */
7042 /* OUTPUT ARGUMENTS : */
7043 /* ------------------- */
7044 /* TABVAL(NDIM): The resulting point (or table of values) */
7046 /* COMMONS USED : */
7047 /* ---------------- */
7049 /* REFERENCES CALLED : */
7050 /* ----------------------- */
7052 /* DESCRIPTION/NOTES/LIMITATIONS : */
7053 /* ----------------------------------- */
7056 /* ***********************************************************************
7059 /* Parameter adjustments */
7061 courbe_dim1 = *ncofmx;
7062 courbe_offset = courbe_dim1 + 1;
7063 courbe -= courbe_offset;
7068 for (nd = 1; nd <= i__1; ++nd) {
7074 for (nd = 1; nd <= i__1; ++nd) {
7075 fu = courbe[*ndeg + nd * courbe_dim1];
7076 for (i__ = *ndeg - 1; i__ >= 1; --i__) {
7077 fu = fu * *tparam + courbe[i__ + nd * courbe_dim1];
7087 //=======================================================================
7088 //function : mmpojac_
7090 //=======================================================================
7091 int mmpojac_(doublereal *tparam,
7101 /* Initialized data */
7105 /* System generated locals */
7106 integer valjac_dim1, i__1, i__2;
7108 /* Local variables */
7109 doublereal cofa, cofb, denom, tnorm[100];
7110 integer ii, jj, kk1, kk2;
7111 doublereal aux1, aux2;
7114 /* ***********************************************************************
7119 /* Positioning on Jacobi polynoms and their derivatives */
7120 /* successive by a recurrent algorithm */
7124 /* RESERVE, POSITIONING, JACOBI */
7126 /* INPUT ARGUMENTS : */
7127 /* -------------------- */
7128 /* TPARAM : Parameter for which positioning is done. */
7129 /* IORDRE : Order of hermit-?? (-1,0,1, or 2) */
7130 /* NCOEFF : Number of coeeficients of polynoms (Nb of value to */
7132 /* NDERIV : Number of derivative to calculate (0<= N <=3) */
7133 /* 0 -> Position simple on jacobi functions */
7134 /* N -> Position on jacobi functions and their */
7135 /* derivatives of order 1 to N. */
7137 /* OUTPUT ARGUMENTS : */
7138 /* --------------------- */
7139 /* VALJAC (NCOEFF, 0:NDERIV) : the calculated values */
7141 /* d vj(t) = VALJAC(J, I) */
7145 /* IERCOD : Error Code */
7147 /* 1 : Incoherence of input arguments */
7149 /* COMMONS USED : */
7150 /* ------------------ */
7153 /* REFERENCES CALLED : */
7154 /* --------------------- */
7157 /* DESCRIPTION/NOTES/LIMITATIONS : */
7158 /* ----------------------------------- */
7161 /* ***********************************************************************
7164 /* ***********************************************************************
7168 /* static varaibles */
7172 /* Parameter adjustments */
7173 valjac_dim1 = *ncoeff;
7178 /* ***********************************************************************
7180 /* INITIALISATIONS */
7181 /* ***********************************************************************
7186 /* ***********************************************************************
7189 /* ***********************************************************************
7195 if (*ncoeff > 100) {
7199 /* --- Calculation of norms */
7201 /* IF (NCOEFF.GT.NBCOF) THEN */
7203 for (ii = 1; ii <= i__1; ++ii) {
7207 for (jj = 1; jj <= i__2; ++jj) {
7208 aux2 = aux2 * (doublereal) (kk1 + *iordre + jj) / (doublereal) (
7211 i__2 = (*iordre << 1) + 1;
7212 tnorm[ii - 1] = sqrt(aux2 * (kk1 * 2. + (*iordre << 1) + 1) / pow__ii(&
7220 /* --- Trivial Positions ----- */
7223 aux1 = (doublereal) (*iordre + 1);
7224 valjac[2] = aux1 * *tparam;
7227 valjac[valjac_dim1 + 1] = 0.;
7228 valjac[valjac_dim1 + 2] = aux1;
7231 valjac[(valjac_dim1 << 1) + 1] = 0.;
7232 valjac[(valjac_dim1 << 1) + 2] = 0.;
7235 valjac[valjac_dim1 * 3 + 1] = 0.;
7236 valjac[valjac_dim1 * 3 + 2] = 0.;
7241 /* --- Positioning by recurrence */
7244 for (ii = 3; ii <= i__1; ++ii) {
7248 aux1 = (doublereal) (*iordre + kk2);
7250 cofa = aux2 * (aux2 + 1) * (aux2 + 2);
7251 cofb = (aux2 + 2) * -2. * aux1 * aux1;
7252 denom = kk1 * 2. * (kk2 + (*iordre << 1) + 1) * aux2;
7256 valjac[ii] = (cofa * *tparam * valjac[kk1] + cofb * valjac[kk2]) *
7260 valjac[ii + valjac_dim1] = (cofa * *tparam * valjac[kk1 +
7261 valjac_dim1] + cofa * valjac[kk1] + cofb * valjac[kk2 +
7262 valjac_dim1]) * denom;
7265 valjac[ii + (valjac_dim1 << 1)] = (cofa * *tparam * valjac[
7266 kk1 + (valjac_dim1 << 1)] + cofa * 2 * valjac[kk1 +
7267 valjac_dim1] + cofb * valjac[kk2 + (valjac_dim1 << 1)]
7272 valjac[ii + valjac_dim1 * 3] = (cofa * *tparam * valjac[kk1 +
7273 valjac_dim1 * 3] + cofa * 3 * valjac[kk1 + (
7274 valjac_dim1 << 1)] + cofb * valjac[kk2 + valjac_dim1 *
7280 /* ---> Normalization */
7283 for (ii = 1; ii <= i__1; ++ii) {
7285 for (jj = 0; jj <= i__2; ++jj) {
7286 valjac[ii + jj * valjac_dim1] = tnorm[ii - 1] * valjac[ii + jj *
7293 /* ***********************************************************************
7295 /* PROCESSING OF ERRORS */
7296 /* ***********************************************************************
7304 /* ***********************************************************************
7306 /* RETURN CALLING PROGRAM */
7307 /* ***********************************************************************
7313 AdvApp2Var_SysBase::maermsg_("MMPOJAC", iercod, 7L);
7318 //=======================================================================
7319 //function : AdvApp2Var_MathBase::mmposui_
7321 //=======================================================================
7322 int AdvApp2Var_MathBase::mmposui_(integer *dimmat,
7329 /* System generated locals */
7332 /* Local variables */
7334 integer imin, jmin, i__, j, k;
7337 /* ***********************************************************************
7342 /* FILL THE TABLE OF POSITIONING POSUIV WHICH ALLOWS TO */
7343 /* PARSE BY COLUMN THE INFERIOR TRIANGULAR PART OF THE */
7344 /* MATRIX IN FORM OF PROFILE */
7349 /* RESERVE, MATRIX, PROFILE */
7351 /* INPUT ARGUMENTS : */
7352 /* -------------------- */
7354 /* NISTOC: NUMBER OF COEFFICIENTS IN THE PROFILE */
7355 /* DIMMAT: NUMBER OF LINE OF THE SYMMETRIC SQUARE MATRIX */
7356 /* APOSIT: TABLE OF POSITIONING OF STORAGE TERMS */
7357 /* APOSIT(1,I) CONTAINS THE NUMBER OF TERMES-1 ON LINE */
7358 /* I IN THE PROFILE OF THE MATRIX */
7359 /* APOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF DIAGONAL TERM */
7363 /* OUTPUT ARGUMENTS : */
7364 /* --------------------- */
7365 /* POSUIV: POSUIV(K) (WHERE K IS THE INDEX OF STORAGE OF MAT(I,J)) */
7366 /* CONTAINS THE SMALLEST NUMBER IMIN>I OF THE LINE THAT */
7367 /* POSSESSES A TERM MAT(IMIN,J) THAT IS IN THE PROFILE. */
7368 /* IF THERE IS NO TERM MAT(IMIN,J) IN THE PROFILE THEN POSUIV(K)=-1 */
7371 /* COMMONS USED : */
7372 /* ------------------ */
7375 /* REFERENCES CALLED : */
7376 /* --------------------- */
7379 /* DESCRIPTION/NOTES/LIMITATIONS : */
7380 /* ----------------------------------- */
7383 /* ***********************************************************************
7386 /* ***********************************************************************
7391 /* ***********************************************************************
7393 /* INITIALIZATIONS */
7394 /* ***********************************************************************
7397 /* Parameter adjustments */
7402 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
7404 AdvApp2Var_SysBase::mgenmsg_("MMPOSUI", 7L);
7409 /* ***********************************************************************
7412 /* ***********************************************************************
7418 for (i__ = 1; i__ <= i__1; ++i__) {
7419 jmin = i__ - aposit[(i__ << 1) + 1];
7421 for (j = jmin; j <= i__2; ++j) {
7424 while(! trouve && imin <= *dimmat) {
7425 if (imin - aposit[(imin << 1) + 1] <= j) {
7431 k = aposit[(i__ << 1) + 2] - i__ + j;
7446 /* ***********************************************************************
7448 /* ERROR PROCESSING */
7449 /* ***********************************************************************
7455 /* ***********************************************************************
7457 /* RETURN CALLING PROGRAM */
7458 /* ***********************************************************************
7463 /* ___ DESALLOCATION, ... */
7465 AdvApp2Var_SysBase::maermsg_("MMPOSUI", iercod, 7L);
7467 AdvApp2Var_SysBase::mgsomsg_("MMPOSUI", 7L);
7472 //=======================================================================
7473 //function : AdvApp2Var_MathBase::mmresol_
7475 //=======================================================================
7476 int AdvApp2Var_MathBase::mmresol_(integer *hdimen,
7494 integer c__100 = 100;
7496 /* System generated locals */
7499 /* Local variables */
7501 doublereal* mcho = 0;
7502 integer jmin, jmax, i__, j, k, l;
7503 intptr_t iofv1, iofv2, iofv3, iofv4;
7504 doublereal *v1 = 0, *v2 = 0, *v3 = 0, *v4 = 0;
7505 integer deblig, dimhch;
7506 doublereal* hchole = 0;
7507 intptr_t iofmch, iofmam, iofhch;
7508 doublereal* matsym = 0;
7514 /* ***********************************************************************
7519 /* SOLUTION OF THE SYSTEM */
7526 /* RESERVE, SOLUTION, SYSTEM, LAGRANGIAN */
7528 /* INPUT ARGUMENTS : */
7529 /* -------------------- */
7530 /* HDIMEN: NOMBER OF LINE (OR COLUMN) OF THE HESSIAN MATRIX */
7531 /* GDIMEN: NOMBER OF LINE OF THE MATRIX OF CONSTRAINTS */
7532 /* HNSTOC: NOMBErS OF TERMS IN THE PROFILE OF HESSIAN MATRIX
7534 /* GNSTOC: NOMBERS OF TERMS IN THE PROFILE OF THE MATRIX OF CONSTRAINTS */
7535 /* MNSTOC: NOMBERS OF TERMS IN THE PROFILE OF THE MATRIX M= G H t(G) */
7536 /* where H IS THE HESSIAN MATRIX AND G IS THE MATRIX OF CONSTRAINTS */
7537 /* MATSYH: TRIANGULAR INFERIOR PART OF THE HESSIAN MATRIX */
7538 /* IN FORM OF PROFILE */
7539 /* MATSYG: MATRIX OF CONSTRAINTS IN FORM OF PROFILE */
7540 /* VECSYH: VECTOR OF THE SECOND MEMBER ASSOCIATED TO MATSYH */
7541 /* VECSYG: VECTOR OF THE SECOND MEMBER ASSOCIATED TO MATSYG */
7542 /* HPOSIT: TABLE OF POSITIONING OF THE HESSIAN MATRIX */
7543 /* HPOSIT(1,I) CONTAINS THE NUMBER OF TERMS -1 */
7544 /* WHICH ARE IN THE PROFILE AT LINE I */
7545 /* HPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF TERM */
7546 /* DIAGONAL OF THE MATRIX AT LINE I */
7547 /* HPOSUI: TABLE ALLOWING TO PARSE THE HESSIAN MATRIX BY COLUMN */
7548 /* IN FORM OF PROFILE */
7549 /* HPOSUI(K) CONTAINS THE NUMBER OF LINE IMIN FOLLOWING THE CURRENT LINE*/
7550 /* I WHERE H(I,J)=MATSYH(K) AS IT EXISTS IN THE */
7551 /* SAME COLUMN J A TERM IN THE PROFILE OF LINE IMIN */
7552 /* IF SUCH TERM DOES NOT EXIST IMIN=-1 */
7553 /* GPOSIT: TABLE OF POSITIONING OF THE MATRIX OF CONSTRAINTS */
7554 /* GPOSIT(1,I) CONTAINS THE NUMBER OF TERMS OF LINE I */
7555 /* WHICH ARE IN THE PROFILE */
7556 /* GPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF THE LAST TERM */
7557 /* OF LINE I WHICH IS IN THE PROFILE */
7558 /* GPOSIT(3,I) CONTAINS THE NUMBER OF COLUMN CORRESPONDING */
7559 /* TO THE FIRST TERM OF LINE I WHICH IS IN THE PROFILE */
7560 /* MMPOSUI, MPOSIT: SAME STRUCTURE AS HPOSUI, BUT FOR MATRIX */
7564 /* OUTPUT ARGUMENTS : */
7565 /* --------------------- */
7566 /* VECSOL: VECTOR SOLUTION V OF THE SYSTEM */
7567 /* IERCOD: ERROR CODE */
7569 /* COMMONS USED : */
7570 /* ------------------ */
7573 /* REFERENCES CALLED : */
7574 /* --------------------- */
7577 /* DESCRIPTION/NOTES/LIMITATIONS : */
7578 /* ----------------------------------- */
7580 /* ***********************************************************************
7583 /* ***********************************************************************
7586 /* ***********************************************************************
7588 /* INITIALISATIONS */
7589 /* ***********************************************************************
7592 /* Parameter adjustments */
7605 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
7607 AdvApp2Var_SysBase::mgenmsg_("MMRESOL", 7L);
7618 /* ***********************************************************************
7621 /* ***********************************************************************
7624 /* Dynamic allocation */
7625 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
7626 anAdvApp2Var_SysBase.macrar8_(hdimen, &c__100, v1, &iofv1, &ier);
7630 dimhch = hposit[(*hdimen << 1) + 2];
7631 anAdvApp2Var_SysBase.macrar8_(&dimhch, &c__100, hchole, &iofhch, &ier);
7636 /* solution of system 1 H V1 = b */
7637 /* where H=MATSYH and b=VECSYH */
7639 mmchole_(hnstoc, hdimen, &matsyh[1], &hposit[3], &hposui[1], &hchole[
7644 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1], &vecsyh[
7645 1], &v1[iofv1], &ier);
7650 /* Case when there are constraints */
7654 /* Calculate the vector of the second member V2=G H(-1) b -c = G v1-c */
7655 /* of system of unknown Lagrangian vector MULTIP */
7656 /* where G=MATSYG */
7659 anAdvApp2Var_SysBase.macrar8_(gdimen, &c__100, v2, &iofv2, &ier);
7663 anAdvApp2Var_SysBase.macrar8_(hdimen, &c__100, v3, &iofv3, &ier);
7667 anAdvApp2Var_SysBase.macrar8_(gdimen, &c__100, v4, &iofv4, &ier);
7671 anAdvApp2Var_SysBase.macrar8_(mnstoc, &c__100, matsym, &iofmam, &ier);
7677 mmatvec_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v1[iofv1], &
7678 deblig, &v2[iofv2], &ier);
7683 for (i__ = 1; i__ <= i__1; ++i__) {
7684 v2[i__ + iofv2 - 1] -= vecsyg[i__];
7687 /* Calculate the matrix M= G H(-1) t(G) */
7688 /* RESOL DU SYST 2 : H qi = gi */
7689 /* where is a vector column of t(G) */
7691 /* then calculate G qi */
7692 /* then construct M in form of profile */
7697 for (i__ = 1; i__ <= i__1; ++i__) {
7698 AdvApp2Var_SysBase::mvriraz_(hdimen, &v1[iofv1]);
7699 AdvApp2Var_SysBase::mvriraz_(hdimen, &v3[iofv3]);
7700 AdvApp2Var_SysBase::mvriraz_(gdimen, &v4[iofv4]);
7701 jmin = gposit[i__ * 3 + 3];
7702 jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
7703 aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
7705 for (j = jmin; j <= i__2; ++j) {
7707 v1[j + iofv1 - 1] = matsyg[k];
7709 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1],
7710 &v1[iofv1], &v3[iofv3], &ier);
7716 mmatvec_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v3[
7717 iofv3], &deblig, &v4[iofv4], &ier);
7722 k = mposit[(i__ << 1) + 2];
7723 matsym[k + iofmam - 1] = v4[i__ + iofv4 - 1];
7724 while(mmposui[k] > 0) {
7726 k = mposit[(l << 1) + 2] - l + i__;
7727 matsym[k + iofmam - 1] = v4[l + iofv4 - 1];
7732 /* SOLVE SYST 3 M L = V2 */
7736 AdvApp2Var_SysBase::mvriraz_(gdimen, &v4[iofv4]);
7737 anAdvApp2Var_SysBase.macrar8_(mnstoc, &c__100, mcho, &iofmch, &ier);
7741 mmchole_(mnstoc, gdimen, &matsym[iofmam], &mposit[3], &mmposui[1], &
7742 mcho[iofmch], &ier);
7746 mmrslss_(mnstoc, gdimen, &mcho[iofmch], &mposit[3], &mmposui[1], &v2[
7747 iofv2], &v4[iofv4], &ier);
7753 /* CALCULATE THE VECTOR OF THE SECOND MEMBER OF THE SYSTEM Hx = b - t(G) L
7757 AdvApp2Var_SysBase::mvriraz_(hdimen, &v1[iofv1]);
7758 mmtmave_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v4[iofv4], &
7764 for (i__ = 1; i__ <= i__1; ++i__) {
7765 v1[i__ + iofv1 - 1] = vecsyh[i__] - v1[i__ + iofv1 - 1];
7768 /* RESOL SYST 4 Hx = b - t(G) L */
7771 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1], &v1[
7772 iofv1], &vecsol[1], &ier);
7778 for (i__ = 1; i__ <= i__1; ++i__) {
7779 vecsol[i__] = v1[i__ + iofv1 - 1];
7785 /* ***********************************************************************
7787 /* PROCESSING OF ERRORS */
7788 /* ***********************************************************************
7797 AdvApp2Var_SysBase::mswrdbg_("MMRESOL : PROBLEM WITH DIMMAT", 30L);
7800 /* ***********************************************************************
7802 /* RETURN CALLING PROGRAM */
7803 /* ***********************************************************************
7808 /* ___ DESALLOCATION, ... */
7809 anAdvApp2Var_SysBase.macrdr8_(hdimen, &c__100, v1, &iofv1, &ier);
7810 if (*iercod == 0 && ier > 0) {
7813 anAdvApp2Var_SysBase.macrdr8_(&dimhch, &c__100, hchole, &iofhch, &ier);
7814 if (*iercod == 0 && ier > 0) {
7817 anAdvApp2Var_SysBase.macrdr8_(gdimen, &c__100, v2, &iofv2, &ier);
7818 if (*iercod == 0 && ier > 0) {
7821 anAdvApp2Var_SysBase.macrdr8_(hdimen, &c__100, v3, &iofv3, &ier);
7822 if (*iercod == 0 && ier > 0) {
7825 anAdvApp2Var_SysBase.macrdr8_(gdimen, &c__100, v4, &iofv4, &ier);
7826 if (*iercod == 0 && ier > 0) {
7829 anAdvApp2Var_SysBase.macrdr8_(mnstoc, &c__100, matsym, &iofmam, &ier);
7830 if (*iercod == 0 && ier > 0) {
7833 anAdvApp2Var_SysBase.macrdr8_(mnstoc, &c__100, mcho, &iofmch, &ier);
7834 if (*iercod == 0 && ier > 0) {
7838 AdvApp2Var_SysBase::maermsg_("MMRESOL", iercod, 7L);
7840 AdvApp2Var_SysBase::mgsomsg_("MMRESOL", 7L);
7845 //=======================================================================
7846 //function : mmrslss_
7848 //=======================================================================
7849 int mmrslss_(integer *,//mxcoef,
7854 doublereal *mscnmbr,
7858 /* System generated locals */
7861 /* Local variables */
7865 integer pointe, ptcour;
7867 /* ***********************************************************************
7872 /* Solves linear system SS x = b where S is a */
7873 /* triangular lower matrix given in form of profile */
7877 /* RESERVE, MATRICE_PROFILE, RESOLUTION, CHOLESKI */
7879 /* INPUT ARGUMENTS : */
7880 /* -------------------- */
7881 /* MXCOEF : Maximum number of non-null coefficient in the matrix */
7882 /* DIMENS : Dimension of the matrix */
7883 /* SMATRI(MXCOEF) : Values of coefficients of the matrix */
7884 /* SPOSIT(2,DIMENS): */
7885 /* SPOSIT(1,*) : Distance diagonal-extremity of the line */
7886 /* SPOSIT(2,*) : Position of diagonal terms in AMATRI */
7887 /* POSUIV(MXCOEF): first line inferior not out of profile */
7888 /* MSCNMBR(DIMENS): Vector second member of the equation */
7890 /* OUTPUT ARGUMENTS : */
7891 /* --------------------- */
7892 /* SOLUTI(NDIMEN) : Result vector */
7893 /* IERCOD : Error code 0 : ok */
7895 /* COMMONS USED : */
7896 /* ------------------ */
7899 /* REFERENCES CALLED : */
7900 /* --------------------- */
7903 /* DESCRIPTION/NOTES/LIMITATIONS : */
7904 /* ----------------------------------- */
7906 /* SS is the decomposition of choleski of a symmetric matrix */
7907 /* defined postive, that can result from routine MMCHOLE. */
7909 /* For a full matrix it is possible to use MRSLMSC */
7911 /* LEVEL OF DEBUG = 4 */
7913 /* ***********************************************************************
7916 /* ***********************************************************************
7921 /* ***********************************************************************
7923 /* INITIALISATIONS */
7924 /* ***********************************************************************
7927 /* Parameter adjustments */
7935 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
7937 AdvApp2Var_SysBase::mgenmsg_("MMRSLSS", 7L);
7941 /* ***********************************************************************
7944 /* ***********************************************************************
7947 /* ----- Solution of Sw = b */
7950 for (i__ = 1; i__ <= i__1; ++i__) {
7952 pointe = sposit[(i__ << 1) + 2];
7955 for (j = i__ - sposit[(i__ << 1) + 1]; j <= i__2; ++j) {
7956 somme += smatri[pointe - (i__ - j)] * soluti[j];
7959 soluti[i__] = (mscnmbr[i__] - somme) / smatri[pointe];
7962 /* ----- Solution of S u = w */
7964 for (i__ = *dimens; i__ >= 1; --i__) {
7966 pointe = sposit[(i__ << 1) + 2];
7970 ptcour = sposit[(j << 1) + 2] - (j - i__);
7971 somme += smatri[ptcour] * soluti[j];
7975 soluti[i__] = (soluti[i__] - somme) / smatri[pointe];
7980 /* ***********************************************************************
7982 /* ERROR PROCESSING */
7983 /* ***********************************************************************
7987 /* ***********************************************************************
7989 /* RETURN PROGRAM CALLING */
7990 /* ***********************************************************************
7995 AdvApp2Var_SysBase::maermsg_("MMRSLSS", iercod, 7L);
7997 AdvApp2Var_SysBase::mgsomsg_("MMRSLSS", 7L);
8002 //=======================================================================
8003 //function : mmrslw_
8005 //=======================================================================
8006 int mmrslw_(integer *normax,
8014 /* System generated locals */
8015 integer abmatr_dim1, abmatr_offset, xmatri_dim1, xmatri_offset, i__1,
8019 /* Local variables */
8026 /* **********************************************************************
8031 /* Solution of a linear system A.x = B of N equations to N */
8032 /* unknown by Gauss method (partial pivot) or : */
8033 /* A is matrix NORDRE * NORDRE, */
8034 /* B is matrix NORDRE (lines) * NDIMEN (columns), */
8035 /* x is matrix NORDRE (lines) * NDIMEN (columns). */
8036 /* In this program, A and B are stored in matrix ABMATR */
8037 /* the lines and columns which of were inverted. ABMATR(k,j) is */
8038 /* term A(j,k) if k <= NORDRE, B(j,k-NORDRE) otherwise (see example). */
8042 /* TOUS, MATH_ACCES::EQUATION&, MATRICE&, RESOLUTION, GAUSS, &SOLUTION */
8044 /* INPUT ARGUMENTS : */
8045 /* ------------------ */
8046 /* NORMAX : Max size of the first index of XMATRI. This argument */
8047 /* serves only for the declaration of dimension of XMATRI and should be */
8048 /* above or equal to NORDRE. */
8049 /* NORDRE : Order of the matrix i.e. number of equations and */
8050 /* unknown quantities of the linear system to be solved. */
8051 /* NDIMEN : Number of the second member. */
8052 /* EPSPIV : Minimal value of a pivot. If during the calculation */
8053 /* the absolute value of the pivot is below EPSPIV, the */
8054 /* system of equations is declared singular. EPSPIV should */
8055 /* be a "small" real. */
8057 /* ABMATR(NORDRE+NDIMEN,NORDRE) : Auxiliary matrix containing */
8058 /* matrix A and matrix B. */
8060 /* OUTPUT ARGUMENTS : */
8061 /* ------------------- */
8062 /* XMATRI : Matrix containing NORDRE*NDIMEN solutions. */
8063 /* IERCOD=0 shows that all solutions are calculated. */
8064 /* IERCOD=1 shows that the matrix is of lower rank than NORDRE */
8065 /* (the system is singular). */
8067 /* COMMONS USED : */
8068 /* ---------------- */
8070 /* REFERENCES CALLED : */
8071 /* ----------------------- */
8073 /* DESCRIPTION/NOTES/LIMITATIONS : */
8074 /* ----------------------------------- */
8075 /* ATTENTION : the indices of line and column are inverted */
8076 /* compared to usual indices. */
8078 /* a1*x + b1*y = c1 */
8079 /* a2*x + b2*y = c2 */
8080 /* should be represented by matrix ABMATR : */
8082 /* ABMATR(1,1) = a1 ABMATR(1,2) = a2 */
8083 /* ABMATR(2,1) = b1 ABMATR(2,2) = b2 */
8084 /* ABMATR(3,1) = c1 ABMATR(3,2) = c2 */
8086 /* To solve this system, it is necessary to set : */
8088 /* NORDRE = 2 (there are 2 equations with 2 unknown values), */
8089 /* NDIMEN = 1 (there is only one second member), */
8090 /* any NORMAX can be taken >= NORDRE. */
8092 /* To use this routine, it is recommended to use one of */
8093 /* interfaces : MMRSLWI or MMMRSLWD. */
8095 /* **********************************************************************
8098 /* Name of the routine */
8100 /* INTEGER IBB,MNFNDEB */
8103 /* IF (IBB.GE.2) CALL MGENMSG(NOMPR) */
8104 /* Parameter adjustments */
8105 xmatri_dim1 = *normax;
8106 xmatri_offset = xmatri_dim1 + 1;
8107 xmatri -= xmatri_offset;
8108 abmatr_dim1 = *nordre + *ndimen;
8109 abmatr_offset = abmatr_dim1 + 1;
8110 abmatr -= abmatr_offset;
8115 /* *********************************************************************
8117 /* Triangulation of matrix ABMATR. */
8118 /* *********************************************************************
8122 for (kk = 1; kk <= i__1; ++kk) {
8124 /* ---------- Find max pivot in column KK. ------------
8130 for (jj = kk; jj <= i__2; ++jj) {
8131 akj = (d__1 = abmatr[kk + jj * abmatr_dim1], advapp_abs(d__1));
8142 /* --------- Swapping of line KPIV with line KK. ------
8146 i__2 = *nordre + *ndimen;
8147 for (jj = kk; jj <= i__2; ++jj) {
8148 akj = abmatr[jj + kk * abmatr_dim1];
8149 abmatr[jj + kk * abmatr_dim1] = abmatr[jj + kpiv *
8151 abmatr[jj + kpiv * abmatr_dim1] = akj;
8156 /* ---------- Removal and triangularization. -----------
8159 pivot = -abmatr[kk + kk * abmatr_dim1];
8161 for (ii = kk + 1; ii <= i__2; ++ii) {
8162 akj = abmatr[kk + ii * abmatr_dim1] / pivot;
8163 i__3 = *nordre + *ndimen;
8164 for (jj = kk + 1; jj <= i__3; ++jj) {
8165 abmatr[jj + ii * abmatr_dim1] += akj * abmatr[jj + kk *
8176 /* *********************************************************************
8178 /* Solution of the system of triangular equations. */
8179 /* Matrix ABMATR(NORDRE+JJ,II), contains second members */
8180 /* of the system for 1<=j<=NDIMEN and 1<=i<=NORDRE. */
8181 /* *********************************************************************
8185 /* ---------------- Calculation of solutions by ascending. -----------------
8188 for (kk = *nordre; kk >= 1; --kk) {
8189 pivot = abmatr[kk + kk * abmatr_dim1];
8191 for (ii = 1; ii <= i__1; ++ii) {
8192 akj = abmatr[ii + *nordre + kk * abmatr_dim1];
8194 for (jj = kk + 1; jj <= i__2; ++jj) {
8195 akj -= abmatr[jj + kk * abmatr_dim1] * xmatri[jj + ii *
8199 xmatri[kk + ii * xmatri_dim1] = akj / pivot;
8206 /* ------If the absolute value of a pivot is smaller than -------- */
8207 /* ---------- EPSPIV: return the code of error. ------------
8217 AdvApp2Var_SysBase::maermsg_("MMRSLW ", iercod, 7L);
8219 /* IF (IBB.GE.2) CALL MGSOMSG(NOMPR) */
8223 //=======================================================================
8224 //function : AdvApp2Var_MathBase::mmmrslwd_
8226 //=======================================================================
8227 int AdvApp2Var_MathBase::mmmrslwd_(integer *normax,
8238 /* System generated locals */
8239 integer amat_dim1, amat_offset, bmat_dim1, bmat_offset, xmat_dim1,
8240 xmat_offset, aaux_dim1, aaux_offset, i__1, i__2;
8242 /* Local variables */
8246 /* IMPLICIT DOUBLE PRECISION (A-H,O-Z) */
8247 /* IMPLICIT INTEGER (I-N) */
8250 /* **********************************************************************
8255 /* Solution of a linear system by Gauss method where */
8256 /* the second member is a table of vectors. Method of partial pivot. */
8260 /* ALL, MATH_ACCES :: */
8261 /* SYSTEME&,EQUATION&, RESOLUTION,GAUSS ,&VECTEUR */
8263 /* INPUT ARGUMENTS : */
8264 /* ------------------ */
8265 /* NORMAX : Max. Dimension of AMAT. */
8266 /* NORDRE : Order of the matrix. */
8267 /* NDIM : Number of columns of BMAT and XMAT. */
8268 /* AMAT(NORMAX,NORDRE) : The processed matrix. */
8269 /* BMAT(NORMAX,NDIM) : The matrix of second member. */
8270 /* XMAT(NORMAX,NDIM) : The matrix of solutions. */
8271 /* EPSPIV : Min value of a pivot. */
8273 /* OUTPUT ARGUMENTS : */
8274 /* ------------------- */
8275 /* AAUX(NORDRE+NDIM,NORDRE) : Auxiliary matrix. */
8276 /* XMAT(NORMAX,NDIM) : Matrix of solutions. */
8277 /* IERCOD=0 shows that solutions in XMAT are valid. */
8278 /* IERCOD=1 shows that matrix AMAT is of lower rank than NORDRE. */
8280 /* COMMONS USED : */
8281 /* ---------------- */
8285 /* REFERENCES CALLED : */
8286 /* ---------------------- */
8288 /* MAERMSG MGENMSG MGSOMSG */
8289 /* MMRSLW I*4 MNFNDEB */
8291 /* DESCRIPTION/NOTES/LIMITATIONS : */
8292 /* ----------------------------------- */
8293 /* ATTENTION : lines and columns are located in usual order : */
8294 /* 1st index = index line */
8295 /* 2nd index = index column */
8296 /* Example, the system : */
8297 /* a1*x + b1*y = c1 */
8298 /* a2*x + b2*y = c2 */
8299 /* is represented by matrix AMAT : */
8301 /* AMAT(1,1) = a1 AMAT(2,1) = a2 */
8302 /* AMAT(1,2) = b1 AMAT(2,2) = b2 */
8304 /* The first index is the index of line, the second index */
8305 /* is the index of columns (Compare with MMRSLWI which is faster). */
8308 /* **********************************************************************
8311 /* Name of the routine */
8313 /* Parameter adjustments */
8314 amat_dim1 = *normax;
8315 amat_offset = amat_dim1 + 1;
8316 amat -= amat_offset;
8317 xmat_dim1 = *normax;
8318 xmat_offset = xmat_dim1 + 1;
8319 xmat -= xmat_offset;
8320 aaux_dim1 = *nordre + *ndim;
8321 aaux_offset = aaux_dim1 + 1;
8322 aaux -= aaux_offset;
8323 bmat_dim1 = *normax;
8324 bmat_offset = bmat_dim1 + 1;
8325 bmat -= bmat_offset;
8328 ibb = AdvApp2Var_SysBase::mnfndeb_();
8330 AdvApp2Var_SysBase::mgenmsg_("MMMRSLW", 7L);
8333 /* Initialization of the auxiliary matrix. */
8336 for (i__ = 1; i__ <= i__1; ++i__) {
8338 for (j = 1; j <= i__2; ++j) {
8339 aaux[j + i__ * aaux_dim1] = amat[i__ + j * amat_dim1];
8345 /* Second member. */
8348 for (i__ = 1; i__ <= i__1; ++i__) {
8350 for (j = 1; j <= i__2; ++j) {
8351 aaux[j + *nordre + i__ * aaux_dim1] = bmat[i__ + j * bmat_dim1];
8357 /* Solution of the system of equations. */
8359 mmrslw_(normax, nordre, ndim, epspiv, &aaux[aaux_offset], &xmat[
8360 xmat_offset], iercod);
8364 AdvApp2Var_SysBase::maermsg_("MMMRSLW", iercod, 7L);
8367 AdvApp2Var_SysBase::mgsomsg_("MMMRSLW", 7L);
8372 //=======================================================================
8373 //function : AdvApp2Var_MathBase::mmrtptt_
8375 //=======================================================================
8376 int AdvApp2Var_MathBase::mmrtptt_(integer *ndglgd,
8380 integer ideb, nmod2, nsur2, ilong, ibb;
8383 /* **********************************************************************
8388 /* Extracts from Common LDGRTL the STRICTLY positive roots of the */
8389 /* Legendre polynom of degree NDGLGD, for 2 <= NDGLGD <= 61. */
8393 /* TOUS, AB_SPECIFI::COMMON&, EXTRACTION, &RACINE, &LEGENDRE. */
8395 /* INPUT ARGUMENTS : */
8396 /* ------------------ */
8397 /* NDGLGD : Mathematic degree of Legendre polynom. */
8398 /* This degree should be above or equal to 2 and */
8399 /* below or equal to 61. */
8401 /* OUTPUT ARGUMENTS : */
8402 /* ------------------- */
8403 /* RTLEGD : The table of strictly positive roots of */
8404 /* Legendre polynom of degree NDGLGD. */
8406 /* COMMONS USED : */
8407 /* ---------------- */
8409 /* REFERENCES CALLED : */
8410 /* ----------------------- */
8412 /* DESCRIPTION/NOTES/LIMITATIONS : */
8413 /* ----------------------------------- */
8414 /* ATTENTION: the condition on NDEGRE ( 2 <= NDEGRE <= 61) is not */
8415 /* tested. The caller should make the test. */
8418 /* **********************************************************************
8420 /* Nome of the routine */
8423 /* Common MLGDRTL: */
8424 /* This common includes POSITIVE roots of Legendre polynoms */
8425 /* AND the weight of Gauss quadrature formulas on all */
8426 /* POSITIVE roots of Legendre polynoms. */
8429 /* ***********************************************************************
8434 /* The common of Legendre roots. */
8440 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
8441 /* ----------------------------------- */
8444 /* ***********************************************************************
8450 /* ROOTAB : Table of all rotts of Legendre polynoms */
8451 /* between [0,1]. They are ranked for degrees increasing from 2 to 61. */
8452 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
8453 /* The address is the same. */
8454 /* HI0TAB : Table of Legendre interpolators for root x=0 */
8455 /* the polynoms of UNEVEN degree. */
8456 /* RTLTB0 : Table of Li(uk) where uk are roots of a */
8457 /* Legendre polynom of EVEN degree. */
8458 /* RTLTB1 : Table of Li(uk) where uk are roots of a */
8459 /* Legendre polynom of UNEVEN degree. */
8462 /************************************************************************
8464 /* Parameter adjustments */
8468 ibb = AdvApp2Var_SysBase::mnfndeb_();
8470 AdvApp2Var_SysBase::mgenmsg_("MMRTPTT", 7L);
8476 nsur2 = *ndglgd / 2;
8477 nmod2 = *ndglgd % 2;
8480 ideb = nsur2 * (nsur2 - 1) / 2 + 1;
8481 AdvApp2Var_SysBase::mcrfill_(&ilong,
8482 &mlgdrtl_.rootab[ideb + nmod2 * 465 - 1],
8485 /* ----------------------------- The end --------------------------------
8490 AdvApp2Var_SysBase::mgsomsg_("MMRTPTT", 7L);
8495 //=======================================================================
8496 //function : AdvApp2Var_MathBase::mmsrre2_
8498 //=======================================================================
8499 int AdvApp2Var_MathBase::mmsrre2_(doublereal *tparam,
8507 /* System generated locals */
8510 /* Local variables */
8511 integer ideb, ifin, imil, ibb;
8513 /* ***********************************************************************
8519 /* Find the interval corresponding to a valueb given in */
8520 /* increasing order of real numbers with double precision. */
8524 /* TOUS,MATH_ACCES::TABLEAU&,POINT&,CORRESPONDANCE,&RANG */
8526 /* INPUT ARGUMENTS : */
8527 /* ------------------ */
8529 /* TPARAM : Value to be tested. */
8530 /* NBRVAL : Size of TABLEV */
8531 /* TABLEV : Table of reals. */
8532 /* EPSIL : Epsilon of precision */
8534 /* OUTPUT ARGUMENTS : */
8535 /* ------------------- */
8537 /* NUMINT : Number of the interval (between 1 and NBRVAL-1). */
8538 /* ITYPEN : = 0 TPARAM is inside the interval NUMINT */
8539 /* = 1 : TPARAM corresponds to the lower limit of */
8540 /* the provided interval. */
8541 /* = 2 : TPARAM corresponds to the upper limit of */
8542 /* the provided interval. */
8544 /* IERCOD : Error code. */
8546 /* = 1 : TABLEV does not contain enough elements. */
8547 /* = 2 : TPARAM out of limits of TABLEV. */
8549 /* COMMONS USED : */
8550 /* ---------------- */
8552 /* REFERENCES CALLED : */
8553 /* ------------------- */
8555 /* DESCRIPTION/NOTES/LIMITATIONS : */
8556 /* --------------------------------- */
8557 /* There are NBRVAL values in TABLEV which stands for NBRVAL-1 intervals. */
8558 /* One searches the interval containing TPARAM by */
8559 /* dichotomy. Complexity of the algorithm : Log(n)/Log(2).(RBD). */
8561 /* ***********************************************************************
8565 /* Initialisations */
8567 /* Parameter adjustments */
8571 ibb = AdvApp2Var_SysBase::mnfndeb_();
8573 AdvApp2Var_SysBase::mgenmsg_("MMSRRE2", 7L);
8582 /* TABLEV should contain at least two values */
8589 /* TPARAM should be between extreme limits of TABLEV. */
8591 if (*tparam < tablev[1] || *tparam > tablev[*nbrval]) {
8596 /* ----------------------- SEARCH OF THE INTERVAL --------------------
8601 /* Test end of loop (found). */
8603 if (ideb + 1 == ifin) {
8608 /* Find by dichotomy on increasing values of TABLEV. */
8610 imil = (ideb + ifin) / 2;
8611 if (*tparam >= tablev[ideb] && *tparam <= tablev[imil]) {
8619 /* -------------- TEST IF TPARAM IS NOT A VALUE --------- */
8620 /* ------------------------OF TABLEV UP TO EPSIL ----------------------
8624 if ((d__1 = *tparam - tablev[ideb], advapp_abs(d__1)) < *epsil) {
8628 if ((d__1 = *tparam - tablev[ifin], advapp_abs(d__1)) < *epsil) {
8633 /* --------------------------- THE END ----------------------------------
8638 AdvApp2Var_SysBase::maermsg_("MMSRRE2", iercod, 7L);
8641 AdvApp2Var_SysBase::mgsomsg_("MMSRRE2", 7L);
8646 //=======================================================================
8647 //function : mmtmave_
8649 //=======================================================================
8650 int mmtmave_(integer *nligne,
8660 /* System generated locals */
8663 /* Local variables */
8665 integer imin, imax, i__, j, k;
8670 /* ***********************************************************************
8676 /* CREATES PRODUCT G V */
8677 /* WHERE THE MATRIX IS IN FORM OF PROFILE */
8681 /* RESERVE, PRODUCT, MATRIX, PROFILE, VECTOR */
8683 /* INPUT ARGUMENTS : */
8684 /* -------------------- */
8685 /* NLIGNE : NUMBER OF LINE OF THE MATRIX */
8686 /* NCOLON : NOMBER OF COLUMN OF THE MATRIX */
8687 /* GPOSIT: TABLE OF POSITIONING OF TERMS OF STORAGE */
8688 /* GPOSIT(1,I) CONTAINS THE NUMBER of TERMS-1 ON LINE */
8689 /* I IN THE PROFILE OF THE MATRIX */
8690 /* GPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF THE DIAGONAL TERM*/
8692 /* GPOSIT(3,I) CONTAINS THE INDEX COLUMN OF THE FIRST TERM OF */
8693 /* PROFILE OF LINE I */
8694 /* GNSTOC : NOMBER OF TERM IN THE PROFILE OF GMATRI */
8695 /* GMATRI : MATRIX OF CONSTRAINTS IN FORM OF PROFILE */
8696 /* VECIN : INPUT VECTOR */
8698 /* OUTPUT ARGUMENTS : */
8699 /* --------------------- */
8700 /* VECOUT : VECTOR PRODUCT */
8701 /* IERCOD : ERROR CODE */
8704 /* COMMONS USED : */
8705 /* ------------------ */
8708 /* REFERENCES CALLED : */
8709 /* --------------------- */
8712 /* DESCRIPTION/NOTES/LIMITATIONS : */
8713 /* ----------------------------------- */
8715 /* ***********************************************************************
8718 /* ***********************************************************************
8723 /* ***********************************************************************
8725 /* INITIALISATIONS */
8726 /* ***********************************************************************
8729 /* Parameter adjustments */
8736 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
8738 AdvApp2Var_SysBase::mgenmsg_("MMTMAVE", 7L);
8742 /* ***********************************************************************
8745 /* ***********************************************************************
8751 for (i__ = 1; i__ <= i__1; ++i__) {
8754 for (j = 1; j <= i__2; ++j) {
8755 imin = gposit[j * 3 + 3];
8756 imax = gposit[j * 3 + 1] + gposit[j * 3 + 3] - 1;
8757 aux = gposit[j * 3 + 2] - gposit[j * 3 + 1] - imin + 1;
8758 if (imin <= i__ && i__ <= imax) {
8760 somme += gmatri[k] * vecin[j];
8763 vecout[i__] = somme;
8772 /* ***********************************************************************
8774 /* ERROR PROCESSING */
8775 /* ***********************************************************************
8779 /* ***********************************************************************
8781 /* RETURN CALLING PROGRAM */
8782 /* ***********************************************************************
8787 /* ___ DESALLOCATION, ... */
8789 AdvApp2Var_SysBase::maermsg_("MMTMAVE", iercod, 7L);
8791 AdvApp2Var_SysBase::mgsomsg_("MMTMAVE", 7L);
8796 //=======================================================================
8797 //function : mmtrpj0_
8799 //=======================================================================
8800 int mmtrpj0_(integer *ncofmx,
8810 /* System generated locals */
8811 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
8814 /* Local variables */
8816 doublereal bidon, error;
8820 /* ***********************************************************************
8825 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
8826 /* Legendre with a given precision. */
8830 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
8832 /* INPUT ARGUMENTS : */
8833 /* ------------------ */
8834 /* NCOFMX : Max Nb of coeff. of the curve (dimensioning). */
8835 /* NDIMEN : Dimension of the space. */
8836 /* NCOEFF : Degree +1 of the polynom. */
8837 /* EPSI3D : Precision required for the approximation. */
8838 /* CRVLGD : The curve the degree which of it is required to lower. */
8840 /* OUTPUT ARGUMENTS : */
8841 /* ------------------- */
8842 /* EPSTRC : Precision of the approximation. */
8843 /* NCFNEW : Degree +1 of the resulting polynom. */
8845 /* COMMONS USED : */
8846 /* ---------------- */
8848 /* REFERENCES CALLED : */
8849 /* ----------------------- */
8851 /* DESCRIPTION/NOTES/LIMITATIONS : */
8852 /* ----------------------------------- */
8854 /* ***********************************************************************
8858 /* ------- Minimum degree that can be attained : Stop at 1 (RBD) ---------
8861 /* Parameter adjustments */
8863 crvlgd_dim1 = *ncofmx;
8864 crvlgd_offset = crvlgd_dim1 + 1;
8865 crvlgd -= crvlgd_offset;
8869 /* ------------------- Init for error calculation -----------------------
8872 for (i__ = 1; i__ <= i__1; ++i__) {
8879 /* Cutting of coefficients. */
8882 /* ------ Loop on the series of Legendre :NCOEFF --> 2 (RBD) -----------
8885 for (i__ = *ncoeff; i__ >= i__1; --i__) {
8886 /* Factor of renormalization. */
8887 bidon = ((i__ - 1) * 2. + 1.) / 2.;
8888 bidon = sqrt(bidon);
8890 for (nd = 1; nd <= i__2; ++nd) {
8891 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
8895 /* Cutting is stopped if the norm becomes too great. */
8896 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
8897 if (error > *epsi3d) {
8902 /* --- Max error cumulee when the I-th coeff is removed. */
8909 /* --------------------------------- End --------------------------------
8916 //=======================================================================
8917 //function : mmtrpj2_
8919 //=======================================================================
8920 int mmtrpj2_(integer *ncofmx,
8930 /* Initialized data */
8932 static doublereal xmaxj[57] = { .9682458365518542212948163499456,
8933 .986013297183269340427888048593603,
8934 1.07810420343739860362585159028115,
8935 1.17325804490920057010925920756025,
8936 1.26476561266905634732910520370741,
8937 1.35169950227289626684434056681946,
8938 1.43424378958284137759129885012494,
8939 1.51281316274895465689402798226634,
8940 1.5878364329591908800533936587012,
8941 1.65970112228228167018443636171226,
8942 1.72874345388622461848433443013543,
8943 1.7952515611463877544077632304216,
8944 1.85947199025328260370244491818047,
8945 1.92161634324190018916351663207101,
8946 1.98186713586472025397859895825157,
8947 2.04038269834980146276967984252188,
8948 2.09730119173852573441223706382076,
8949 2.15274387655763462685970799663412,
8950 2.20681777186342079455059961912859,
8951 2.25961782459354604684402726624239,
8952 2.31122868752403808176824020121524,
8953 2.36172618435386566570998793688131,
8954 2.41117852396114589446497298177554,
8955 2.45964731268663657873849811095449,
8956 2.50718840313973523778244737914028,
8957 2.55385260994795361951813645784034,
8958 2.59968631659221867834697883938297,
8959 2.64473199258285846332860663371298,
8960 2.68902863641518586789566216064557,
8961 2.73261215675199397407027673053895,
8962 2.77551570192374483822124304745691,
8963 2.8177699459714315371037628127545,
8964 2.85940333797200948896046563785957,
8965 2.90044232019793636101516293333324,
8966 2.94091151970640874812265419871976,
8967 2.98083391718088702956696303389061,
8968 3.02023099621926980436221568258656,
8969 3.05912287574998661724731962377847,
8970 3.09752842783622025614245706196447,
8971 3.13546538278134559341444834866301,
8972 3.17295042316122606504398054547289,
8973 3.2099992681699613513775259670214,
8974 3.24662674946606137764916854570219,
8975 3.28284687953866689817670991319787,
8976 3.31867291347259485044591136879087,
8977 3.35411740487202127264475726990106,
8978 3.38919225660177218727305224515862,
8979 3.42390876691942143189170489271753,
8980 3.45827767149820230182596660024454,
8981 3.49230918177808483937957161007792,
8982 3.5260130200285724149540352829756,
8983 3.55939845146044235497103883695448,
8984 3.59247431368364585025958062194665,
8985 3.62524904377393592090180712976368,
8986 3.65773070318071087226169680450936,
8987 3.68992700068237648299565823810245,
8988 3.72184531357268220291630708234186 };
8990 /* System generated locals */
8991 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
8994 /* Local variables */
8996 doublereal bidon, error;
8998 doublereal bid, eps1;
9001 /* ***********************************************************************
9006 /* Lower the degree of a curve defined on (-1,1) in the direction of */
9007 /* Legendre with a given precision. */
9011 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
9013 /* INPUT ARGUMENTS : */
9014 /* ------------------ */
9015 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9016 /* NDIMEN : Dimension of the space. */
9017 /* NCOEFF : Degree +1 of the polynom. */
9018 /* EPSI3D : Precision required for the approximation. */
9019 /* CRVLGD : The curve the degree which of will be lowered. */
9021 /* OUTPUT ARGUMENTS : */
9022 /* ------------------- */
9023 /* YCVMAX : Auxiliary table (error max on each dimension).
9025 /* EPSTRC : Precision of the approximation. */
9026 /* NCFNEW : Degree +1 of the resulting polynom. */
9028 /* COMMONS USED : */
9029 /* ---------------- */
9031 /* REFERENCES CALLED : */
9032 /* ----------------------- */
9034 /* DESCRIPTION/NOTES/LIMITATIONS : */
9035 /* ----------------------------------- */
9037 /* ***********************************************************************
9041 /* Parameter adjustments */
9043 crvlgd_dim1 = *ncofmx;
9044 crvlgd_offset = crvlgd_dim1 + 1;
9045 crvlgd -= crvlgd_offset;
9051 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9055 /* Init for calculation of error. */
9057 for (i__ = 1; i__ <= i__1; ++i__) {
9064 /* Cutting of coefficients. */
9067 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9070 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9071 /* Factor of renormalization. */
9072 bidon = xmaxj[i__ - ncut];
9074 for (nd = 1; nd <= i__2; ++nd) {
9075 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
9079 /* One stops to cut if the norm becomes too great. */
9080 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9081 if (error > *epsi3d) {
9086 /* --- Max error cumulated when the I-th coeff is removed. */
9093 /* ------- Cutting of zero coeffs of interpolation (RBD) -------
9097 if (*ncfnew == ia) {
9098 AdvApp2Var_MathBase::mmeps1_(&eps1);
9099 for (i__ = ia; i__ >= 2; --i__) {
9102 for (nd = 1; nd <= i__1; ++nd) {
9103 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
9112 /* --- If all coeffs can be removed, this is a point. */
9116 /* --------------------------------- End --------------------------------
9123 //=======================================================================
9124 //function : mmtrpj4_
9126 //=======================================================================
9127 int mmtrpj4_(integer *ncofmx,
9136 /* Initialized data */
9138 static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
9139 1.05299572648705464724876659688996,
9140 1.0949715351434178709281698645813,
9141 1.15078388379719068145021100764647,
9142 1.2094863084718701596278219811869,
9143 1.26806623151369531323304177532868,
9144 1.32549784426476978866302826176202,
9145 1.38142537365039019558329304432581,
9146 1.43575531950773585146867625840552,
9147 1.48850442653629641402403231015299,
9148 1.53973611681876234549146350844736,
9149 1.58953193485272191557448229046492,
9150 1.63797820416306624705258190017418,
9151 1.68515974143594899185621942934906,
9152 1.73115699602477936547107755854868,
9153 1.77604489805513552087086912113251,
9154 1.81989256661534438347398400420601,
9155 1.86276344480103110090865609776681,
9156 1.90471563564740808542244678597105,
9157 1.94580231994751044968731427898046,
9158 1.98607219357764450634552790950067,
9159 2.02556989246317857340333585562678,
9160 2.06433638992049685189059517340452,
9161 2.10240936014742726236706004607473,
9162 2.13982350649113222745523925190532,
9163 2.17661085564771614285379929798896,
9164 2.21280102016879766322589373557048,
9165 2.2484214321456956597803794333791,
9166 2.28349755104077956674135810027654,
9167 2.31805304852593774867640120860446,
9168 2.35210997297725685169643559615022,
9169 2.38568889602346315560143377261814,
9170 2.41880904328694215730192284109322,
9171 2.45148841120796359750021227795539,
9172 2.48374387161372199992570528025315,
9173 2.5155912654873773953959098501893,
9174 2.54704548720896557684101746505398,
9175 2.57812056037881628390134077704127,
9176 2.60882970619319538196517982945269,
9177 2.63918540521920497868347679257107,
9178 2.66919945330942891495458446613851,
9179 2.69888301230439621709803756505788,
9180 2.72824665609081486737132853370048,
9181 2.75730041251405791603760003778285,
9182 2.78605380158311346185098508516203,
9183 2.81451587035387403267676338931454,
9184 2.84269522483114290814009184272637,
9185 2.87060005919012917988363332454033,
9186 2.89823818258367657739520912946934,
9187 2.92561704377132528239806135133273,
9188 2.95274375377994262301217318010209,
9189 2.97962510678256471794289060402033,
9190 3.00626759936182712291041810228171,
9191 3.03267744830655121818899164295959,
9192 3.05886060707437081434964933864149 };
9194 /* System generated locals */
9195 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
9198 /* Local variables */
9200 doublereal bidon, error;
9202 doublereal bid, eps1;
9206 /* ***********************************************************************
9211 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
9212 /* Legendre with a given precision. */
9216 /* LEGENDRE, POLYGON, TRONCATION, CURVE, SMOOTHING. */
9218 /* INPUT ARGUMENTS : */
9219 /* ------------------ */
9220 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9221 /* NDIMEN : Dimension of the space. */
9222 /* NCOEFF : Degree +1 of the polynom. */
9223 /* EPSI3D : Precision required for the approximation. */
9224 /* CRVLGD : The curve which wishes to lower the degree. */
9226 /* OUTPUT ARGUMENTS : */
9227 /* ------------------- */
9228 /* YCVMAX : Auxiliary table (max error on each dimension).
9230 /* EPSTRC : Precision of the approximation. */
9231 /* NCFNEW : Degree +1 of the resulting polynom. */
9233 /* COMMONS USED : */
9234 /* ---------------- */
9236 /* REFERENCES CALLED : */
9237 /* ----------------------- */
9239 /* DESCRIPTION/NOTES/LIMITATIONS : */
9240 /* ----------------------------------- */
9242 /* ***********************************************************************
9246 /* Parameter adjustments */
9248 crvlgd_dim1 = *ncofmx;
9249 crvlgd_offset = crvlgd_dim1 + 1;
9250 crvlgd -= crvlgd_offset;
9256 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9260 /* Init for error calculation. */
9262 for (i__ = 1; i__ <= i__1; ++i__) {
9269 /* Cutting of coefficients. */
9272 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9275 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9276 /* Factor of renormalization. */
9277 bidon = xmaxj[i__ - ncut];
9279 for (nd = 1; nd <= i__2; ++nd) {
9280 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
9284 /* Stop cutting if the norm becomes too great. */
9285 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9286 if (error > *epsi3d) {
9291 /* -- Error max cumulated when the I-eme coeff is removed. */
9298 /* ------- Cutting of zero coeffs of the pole of interpolation (RBD) -------
9302 if (*ncfnew == ia) {
9303 AdvApp2Var_MathBase::mmeps1_(&eps1);
9304 for (i__ = ia; i__ >= 2; --i__) {
9307 for (nd = 1; nd <= i__1; ++nd) {
9308 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
9317 /* --- If all coeffs can be removed, this is a point. */
9321 /* --------------------------------- End --------------------------------
9328 //=======================================================================
9329 //function : mmtrpj6_
9331 //=======================================================================
9332 int mmtrpj6_(integer *ncofmx,
9342 /* Initialized data */
9344 static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
9345 1.11626917091567929907256116528817,
9346 1.1327140810290884106278510474203,
9347 1.1679452722668028753522098022171,
9348 1.20910611986279066645602153641334,
9349 1.25228283758701572089625983127043,
9350 1.29591971597287895911380446311508,
9351 1.3393138157481884258308028584917,
9352 1.3821288728999671920677617491385,
9353 1.42420414683357356104823573391816,
9354 1.46546895108549501306970087318319,
9355 1.50590085198398789708599726315869,
9356 1.54550385142820987194251585145013,
9357 1.58429644271680300005206185490937,
9358 1.62230484071440103826322971668038,
9359 1.65955905239130512405565733793667,
9360 1.69609056468292429853775667485212,
9361 1.73193098017228915881592458573809,
9362 1.7671112206990325429863426635397,
9363 1.80166107681586964987277458875667,
9364 1.83560897003644959204940535551721,
9365 1.86898184653271388435058371983316,
9366 1.90180515174518670797686768515502,
9367 1.93410285411785808749237200054739,
9368 1.96589749778987993293150856865539,
9369 1.99721027139062501070081653790635,
9370 2.02806108474738744005306947877164,
9371 2.05846864831762572089033752595401,
9372 2.08845055210580131460156962214748,
9373 2.11802334209486194329576724042253,
9374 2.14720259305166593214642386780469,
9375 2.17600297710595096918495785742803,
9376 2.20443832785205516555772788192013,
9377 2.2325216999457379530416998244706,
9378 2.2602654243075083168599953074345,
9379 2.28768115912702794202525264301585,
9380 2.3147799369092684021274946755348,
9381 2.34157220782483457076721300512406,
9382 2.36806787963276257263034969490066,
9383 2.39427635443992520016789041085844,
9384 2.42020656255081863955040620243062,
9385 2.44586699364757383088888037359254,
9386 2.47126572552427660024678584642791,
9387 2.49641045058324178349347438430311,
9388 2.52130850028451113942299097584818,
9389 2.54596686772399937214920135190177,
9390 2.5703922285006754089328998222275,
9391 2.59459096001908861492582631591134,
9392 2.61856915936049852435394597597773,
9393 2.64233265984385295286445444361827,
9394 2.66588704638685848486056711408168,
9395 2.68923766976735295746679957665724,
9396 2.71238965987606292679677228666411 };
9398 /* System generated locals */
9399 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
9402 /* Local variables */
9404 doublereal bidon, error;
9406 doublereal bid, eps1;
9410 /* ***********************************************************************
9415 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
9416 /* Legendre to a given precision. */
9420 /* LEGENDRE,POLYGON,TRUNCATION,CURVE,SMOOTHING. */
9422 /* INPUT ARGUMENTS : */
9423 /* ------------------ */
9424 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9425 /* NDIMEN : Dimension of the space. */
9426 /* NCOEFF : Degree +1 of the polynom. */
9427 /* EPSI3D : Precision required for the approximation. */
9428 /* CRVLGD : The curve the degree which of will be lowered. */
9430 /* OUTPUT ARGUMENTS : */
9431 /* ------------------- */
9432 /* YCVMAX : Auxiliary table (max error on each dimension). */
9433 /* EPSTRC : Precision of the approximation. */
9434 /* NCFNEW : Degree +1 of the resulting polynom. */
9436 /* COMMONS USED : */
9437 /* ---------------- */
9439 /* REFERENCES CALLED : */
9440 /* ----------------------- */
9442 /* DESCRIPTION/NOTES/LIMITATIONS : */
9443 /* ----------------------------------- */
9445 /* ***********************************************************************
9449 /* Parameter adjustments */
9451 crvlgd_dim1 = *ncofmx;
9452 crvlgd_offset = crvlgd_dim1 + 1;
9453 crvlgd -= crvlgd_offset;
9459 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9463 /* Init for error calculation. */
9465 for (i__ = 1; i__ <= i__1; ++i__) {
9472 /* Cutting of coefficients. */
9475 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9478 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9479 /* Factor of renormalization. */
9480 bidon = xmaxj[i__ - ncut];
9482 for (nd = 1; nd <= i__2; ++nd) {
9483 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
9487 /* Stop cutting if the norm becomes too great. */
9488 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9489 if (error > *epsi3d) {
9494 /* --- Max error cumulated when the I-th coeff is removed. */
9501 /* ------- Cutting of zero coeff. of the pole of interpolation (RBD) -------
9505 if (*ncfnew == ia) {
9506 AdvApp2Var_MathBase::mmeps1_(&eps1);
9507 for (i__ = ia; i__ >= 2; --i__) {
9510 for (nd = 1; nd <= i__1; ++nd) {
9511 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
9520 /* --- If all coeffs can be removed, this is a point. */
9524 /* --------------------------------- End --------------------------------
9531 //=======================================================================
9532 //function : AdvApp2Var_MathBase::mmtrpjj_
9534 //=======================================================================
9535 int AdvApp2Var_MathBase::mmtrpjj_(integer *ncofmx,
9545 /* System generated locals */
9546 integer crvlgd_dim1, crvlgd_offset;
9548 /* Local variables */
9552 /* ***********************************************************************
9557 /* Lower the degree of a curve defined on (-1,1) in the direction of */
9558 /* Legendre with a given precision. */
9562 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
9564 /* INPUT ARGUMENTS : */
9565 /* ------------------ */
9566 /* NCOFMX : Max Nb coeff. of the curve (dimensioning). */
9567 /* NDIMEN : Dimension of the space. */
9568 /* NCOEFF : Degree +1 of the polynom. */
9569 /* EPSI3D : Precision required for the approximation. */
9570 /* IORDRE : Order of continuity at the extremities. */
9571 /* CRVLGD : The curve the degree which of should be lowered. */
9573 /* OUTPUT ARGUMENTS : */
9574 /* ------------------- */
9575 /* ERRMAX : Precision of the approximation. */
9576 /* NCFNEW : Degree +1 of the resulting polynom. */
9578 /* COMMONS USED : */
9579 /* ---------------- */
9581 /* REFERENCES CALLED : */
9582 /* ----------------------- */
9584 /* DESCRIPTION/NOTES/LIMITATIONS : */
9585 /* ----------------------------------- */
9587 /* ***********************************************************************
9591 /* Parameter adjustments */
9593 crvlgd_dim1 = *ncofmx;
9594 crvlgd_offset = crvlgd_dim1 + 1;
9595 crvlgd -= crvlgd_offset;
9598 ia = (*iordre + 1) << 1;
9601 mmtrpj0_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9602 ycvmax[1], errmax, ncfnew);
9603 } else if (ia == 2) {
9604 mmtrpj2_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9605 ycvmax[1], errmax, ncfnew);
9606 } else if (ia == 4) {
9607 mmtrpj4_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9608 ycvmax[1], errmax, ncfnew);
9610 mmtrpj6_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9611 ycvmax[1], errmax, ncfnew);
9614 /* ------------------------ End -----------------------------------------
9620 //=======================================================================
9621 //function : AdvApp2Var_MathBase::mmunivt_
9623 //=======================================================================
9624 int AdvApp2Var_MathBase::mmunivt_(integer *ndimen,
9631 doublereal c_b2 = 10.;
9633 /* System generated locals */
9637 /* Local variables */
9638 integer nchif, iunit = 1, izero;
9647 /* ***********************************************************************
9652 /* CALCULATE THE NORMAL VECTOR BASING ON ANY VECTOR */
9653 /* WITH PRECISION GIVEN BY THE USER. */
9657 /* ALL, MATH_ACCES :: */
9658 /* VECTEUR&, NORMALISATION, &VECTEUR */
9660 /* INPUT ARGUMENTS : */
9661 /* ------------------ */
9662 /* NDIMEN : DIMENSION OF THE SPACE */
9663 /* VECTOR : VECTOR TO BE NORMED */
9664 /* EPSILN : EPSILON BELOW WHICH IT IS CONSIDERED THAT THE */
9665 /* NORM OF THE VECTOR IS NULL. IF EPSILN<=0, A DEFAULT VALUE */
9666 /* IS IMPOSED (10.D-17 ON VAX). */
9668 /* OUTPUT ARGUMENTS : */
9669 /* ------------------- */
9670 /* VECNRM : NORMED VECTOR */
9671 /* IERCOD 101 : THE VECTOR IS NULL UP TO EPSILN. */
9674 /* COMMONS USED : */
9675 /* ---------------- */
9677 /* REFERENCES CALLED : */
9678 /* ----------------------- */
9680 /* DESCRIPTION/NOTES/LIMITATIONS : */
9681 /* ----------------------------------- */
9682 /* VECTOR and VECNRM can be identic. */
9684 /* The norm of vector is calculated and each component is divided by */
9685 /* this norm. After this it is checked if all componentes of the */
9686 /* vector except for one cost 0 with machine precision. In */
9687 /* this case the quasi-null components are set to 0.D0. */
9689 /* ***********************************************************************
9693 /* Parameter adjustments */
9700 /* -------- Precision by default : zero machine 10.D-17 on Vax ------
9703 AdvApp2Var_SysBase::maovsr8_(&nchif);
9704 if (*epsiln <= 0.) {
9706 eps0 = AdvApp2Var_MathBase::pow__di(&c_b2, &i__1);
9711 /* ------------------------- Calculation of the norm --------------------
9714 vnorm = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vector[1]);
9715 if (vnorm <= eps0) {
9716 AdvApp2Var_SysBase::mvriraz_(ndimen, &vecnrm[1]);
9721 /* ---------------------- Calculation of the vector norm ---------------
9725 i__1 = (-nchif - 1) / 2;
9726 eps0 = AdvApp2Var_MathBase::pow__di(&c_b2, &i__1);
9728 for (ii = 1; ii <= i__1; ++ii) {
9729 vecnrm[ii] = vector[ii] / vnorm;
9730 if ((d__1 = vecnrm[ii], advapp_abs(d__1)) <= eps0) {
9738 /* ------ Case when all coordinates except for one are almost null ----
9740 /* ------------- then one of coordinates costs 1.D0 or -1.D0 --------
9743 if (izero == *ndimen - 1) {
9744 bid = vecnrm[iunit];
9746 for (ii = 1; ii <= i__1; ++ii) {
9753 vecnrm[iunit] = -1.;
9757 /* -------------------------------- The end -----------------------------
9764 //=======================================================================
9765 //function : AdvApp2Var_MathBase::mmveps3_
9767 //=======================================================================
9768 int AdvApp2Var_MathBase::mmveps3_(doublereal *eps03)
9770 /* Initialized data */
9772 static char nomprg[8+1] = "MMEPS1 ";
9778 /************************************************************************
9783 /* Extraction of EPS1 from COMMON MPRCSN. */
9787 /* MPRCSN,PRECISON,EPS3. */
9789 /* INPUT ARGUMENTS : */
9790 /* ------------------ */
9793 /* OUTPUT ARGUMENTS : */
9794 /* ------------------- */
9795 /* EPS3 : space zero of the denominator (10**-9) */
9796 /* EPS3 should value 10**-15 */
9798 /* COMMONS USED : */
9799 /* ---------------- */
9801 /* REFERENCES CALLED : */
9802 /* ----------------------- */
9804 /* DESCRIPTION/NOTES/LIMITATIONS : */
9805 /* ----------------------------------- */
9808 /* ***********************************************************************
9813 /* ***********************************************************************
9818 /* GIVES TOLERANCES OF NULLITY IN STRIM */
9819 /* AND LIMITS OF ITERATIVE PROCESSES */
9821 /* GENERAL CONTEXT, MODIFIABLE BY THE UTILISER */
9825 /* PARAMETER , TOLERANCE */
9827 /* DESCRIPTION/NOTES/LIMITATIONS : */
9828 /* ----------------------------------- */
9829 /* INITIALISATION : PROFILE , **VIA MPRFTX** AT INPUT IN STRIM*/
9830 /* LOADING OF DEFAULT VALUES OF THE PROFILE IN MPRFTX AT INPUT*/
9831 /* IN STRIM. THEY ARE PRESERVED IN THE LOCAL VARIABLES OF MPRFTX */
9833 /* RESET DEFAULT VALUES : MDFINT */
9834 /* MODIFICATION INTERACTIVE BY THE USER : MDBINT */
9836 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
9837 /* MEPSPB ... EPS3,EPS4 */
9838 /* MEPSLN ... EPS2, NITERM , NITERR */
9839 /* MEPSNR ... EPS2 , NITERM */
9840 /* MITERR ... NITERR */
9843 /* ***********************************************************************
9846 /* NITERM : MAX NB OF ITERATIONS */
9847 /* NITERR : NB OF RAPID ITERATIONS */
9848 /* EPS1 : TOLERANCE OF 3D NULL DISTANCE */
9849 /* EPS2 : TOLERANCE OF ZERO PARAMETRIC DISTANCE */
9850 /* EPS3 : TOLERANCE TO AVOID DIVISION BY 0.. */
9851 /* EPS4 : TOLERANCE ANGULAR */
9855 /* ***********************************************************************
9858 ibb = AdvApp2Var_SysBase::mnfndeb_();
9860 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
9863 *eps03 = mmprcsn_.eps3;
9868 //=======================================================================
9869 //function : AdvApp2Var_MathBase::mmvncol_
9871 //=======================================================================
9872 int AdvApp2Var_MathBase::mmvncol_(integer *ndimen,
9878 /* System generated locals */
9881 /* Local variables */
9884 doublereal vaux1[3], vaux2[3];
9890 /* ***********************************************************************
9895 /* CALCULATE A VECTOR NON-COLINEAR TO A GIVEN NON-NULL VECTOR */
9899 /* PUBLIC, VECTOR, FREE */
9901 /* INPUT ARGUMENTS : */
9902 /* -------------------- */
9903 /* ndimen : dimension of the space */
9904 /* vecin : input vector */
9906 /* OUTPUT ARGUMENTS : */
9907 /* --------------------- */
9909 /* vecout : vector non colinear to vecin */
9911 /* COMMONS USED : */
9912 /* ------------------ */
9915 /* REFERENCES CALLED : */
9916 /* --------------------- */
9919 /* DESCRIPTION/NOTES/LIMITATIONS : */
9920 /* ----------------------------------- */
9922 /* ***********************************************************************
9925 /* ***********************************************************************
9930 /* ***********************************************************************
9932 /* INITIALISATIONS */
9933 /* ***********************************************************************
9936 /* Parameter adjustments */
9941 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
9943 AdvApp2Var_SysBase::mgenmsg_("MMVNCOL", 7L);
9947 /* ***********************************************************************
9950 /* ***********************************************************************
9953 if (*ndimen <= 1 || *ndimen > 3) {
9959 while(d__ <= *ndimen) {
9960 if (vecin[d__] == 0.) {
9965 if (aux == *ndimen) {
9970 for (d__ = 1; d__ <= 3; ++d__) {
9971 vaux1[d__ - 1] = 0.;
9974 for (d__ = 1; d__ <= i__1; ++d__) {
9975 vaux1[d__ - 1] = vecin[d__];
9976 vaux2[d__ - 1] = vecin[d__];
9985 vaux2[d__ - 1] += 1;
9986 valaux = vaux1[1] * vaux2[2] - vaux1[2] * vaux2[1];
9988 valaux = vaux1[2] * vaux2[0] - vaux1[0] * vaux2[2];
9990 valaux = vaux1[0] * vaux2[1] - vaux1[1] * vaux2[0];
10005 for (d__ = 1; d__ <= i__1; ++d__) {
10006 vecout[d__] = vaux2[d__ - 1];
10011 /* ***********************************************************************
10013 /* ERROR PROCESSING */
10014 /* ***********************************************************************
10023 /* ***********************************************************************
10025 /* RETURN CALLING PROGRAM */
10026 /* ***********************************************************************
10032 AdvApp2Var_SysBase::maermsg_("MMVNCOL", iercod, 7L);
10034 AdvApp2Var_SysBase::mgsomsg_("MMVNCOL", 7L);
10039 //=======================================================================
10040 //function : AdvApp2Var_MathBase::mmwprcs_
10042 //=======================================================================
10043 void AdvApp2Var_MathBase::mmwprcs_(doublereal *epsil1,
10044 doublereal *epsil2,
10045 doublereal *epsil3,
10046 doublereal *epsil4,
10053 /* ***********************************************************************
10058 /* ACCESS IN WRITING FOR COMMON MPRCSN */
10064 /* INPUT ARGUMENTS : */
10065 /* -------------------- */
10066 /* EPSIL1 : TOLERANCE OF 3D NULL DISTANCE */
10067 /* EPSIL2 : TOLERANCE OF PARAMETRIC NULL DISTANCE */
10068 /* EPSIL3 : TOLERANCE TO AVOID DIVISION BY 0.. */
10069 /* EPSIL4 : ANGULAR TOLERANCE */
10070 /* NITER1 : MAX NB OF ITERATIONS */
10071 /* NITER2 : NB OF RAPID ITERATIONS */
10073 /* OUTPUT ARGUMENTS : */
10074 /* --------------------- */
10077 /* COMMONS USED : */
10078 /* ------------------ */
10081 /* REFERENCES CALLED : */
10082 /* --------------------- */
10085 /* DESCRIPTION/NOTES/LIMITATIONS : */
10086 /* ----------------------------------- */
10089 /* ***********************************************************************
10092 /* ***********************************************************************
10096 /* ***********************************************************************
10098 /* INITIALIZATIONS */
10099 /* ***********************************************************************
10102 /* ***********************************************************************
10105 /* ***********************************************************************
10108 /* ***********************************************************************
10113 /* GIVES TOLERANCES OF NULLITY IN STRIM */
10114 /* AND LIMITS OF ITERATIVE PROCESSES */
10116 /* GENERAL CONTEXT, MODIFIABLE BY THE UTILISER */
10120 /* PARAMETER , TOLERANCE */
10122 /* DESCRIPTION/NOTES/LIMITATIONS : */
10123 /* ----------------------------------- */
10124 /* INITIALISATION : PROFILE , **VIA MPRFTX** AT INPUT IN STRIM*/
10125 /* LOADING OF DEFAULT VALUES OF THE PROFILE IN MPRFTX AT INPUT*/
10126 /* IN STRIM. THEY ARE PRESERVED IN THE LOCAL VARIABLES OF MPRFTX */
10128 /* RESET DEFAULT VALUES : MDFINT */
10129 /* MODIFICATION INTERACTIVE BY THE USER : MDBINT */
10131 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
10132 /* MEPSPB ... EPS3,EPS4 */
10133 /* MEPSLN ... EPS2, NITERM , NITERR */
10134 /* MEPSNR ... EPS2 , NITERM */
10135 /* MITERR ... NITERR */
10138 /* ***********************************************************************
10141 /* NITERM : MAX NB OF ITERATIONS */
10142 /* NITERR : NB OF RAPID ITERATIONS */
10143 /* EPS1 : TOLERANCE OF 3D NULL DISTANCE */
10144 /* EPS2 : TOLERANCE OF ZERO PARAMETRIC DISTANCE */
10145 /* EPS3 : TOLERANCE TO AVOID DIVISION BY 0.. */
10146 /* EPS4 : TOLERANCE ANGULAR */
10149 /* ***********************************************************************
10151 mmprcsn_.eps1 = *epsil1;
10152 mmprcsn_.eps2 = *epsil2;
10153 mmprcsn_.eps3 = *epsil3;
10154 mmprcsn_.eps4 = *epsil4;
10155 mmprcsn_.niterm = *niter1;
10156 mmprcsn_.niterr = *niter2;
10161 //=======================================================================
10162 //function : AdvApp2Var_MathBase::pow__di
10164 //=======================================================================
10165 doublereal AdvApp2Var_MathBase::pow__di (doublereal *x,
10169 register integer ii ;
10170 doublereal result ;
10173 if ( *n > 0 ) {absolute = *n;}
10174 else {absolute = -*n;}
10175 /* System generated locals */
10176 for(ii = 0 ; ii < absolute ; ii++) {
10180 result = 1.0e0 / result ;
10186 /* **********************************************************************
10191 /* Calculate integer function power not obligatory in the most efficient way ;
10198 /* INPUT ARGUMENTS : */
10199 /* ------------------ */
10200 /* X : argument of X**N */
10203 /* OUTPUT ARGUMENTS : */
10204 /* ------------------- */
10207 /* COMMONS USED : */
10208 /* ---------------- */
10210 /* REFERENCES CALLED : */
10211 /* ----------------------- */
10213 /* DESCRIPTION/NOTES/LIMITATIONS : */
10214 /* ----------------------------------- */
10217 /* ***********************************************************************/
10219 //=======================================================================
10220 //function : pow__ii
10222 //=======================================================================
10223 integer pow__ii(integer *x,
10227 register integer ii ;
10231 if ( *n > 0 ) {absolute = *n;}
10232 else {absolute = -*n;}
10233 /* System generated locals */
10234 for(ii = 0 ; ii < absolute ; ii++) {
10238 result = 1 / result ;
10244 /* **********************************************************************
10246 /* **********************************************************************
10251 /* Calculate integer function power not obligatory in the most efficient way ;
10258 /* INPUT ARGUMENTS : */
10259 /* ------------------ */
10260 /* X : argument of X**N */
10263 /* OUTPUT ARGUMENTS : */
10264 /* ------------------- */
10267 /* COMMONS USED : */
10268 /* ---------------- */
10270 /* REFERENCES CALLED : */
10271 /* ----------------------- */
10273 /* DESCRIPTION/NOTES/LIMITATIONS : */
10274 /* ----------------------------------- */
10277 /* ***********************************************************************/
10279 //=======================================================================
10280 //function : AdvApp2Var_MathBase::msc_
10282 //=======================================================================
10283 doublereal AdvApp2Var_MathBase::msc_(integer *ndimen,
10284 doublereal *vecte1,
10285 doublereal *vecte2)
10288 /* System generated locals */
10290 doublereal ret_val;
10292 /* Local variables */
10298 /************************************************************************
10303 /* Calculate the scalar product of 2 vectors in the space */
10304 /* of dimension NDIMEN. */
10308 /* PRODUCT MSCALAIRE. */
10310 /* INPUT ARGUMENTS : */
10311 /* ------------------ */
10312 /* NDIMEN : Dimension of the space. */
10313 /* VECTE1,VECTE2: Vectors. */
10315 /* OUTPUT ARGUMENTS : */
10316 /* ------------------- */
10318 /* COMMONS USED : */
10319 /* ---------------- */
10321 /* REFERENCES CALLED : */
10322 /* ----------------------- */
10324 /* DESCRIPTION/NOTES/LIMITATIONS : */
10325 /* ----------------------------------- */
10328 /* ***********************************************************************
10332 /* PRODUIT MSCALAIRE */
10333 /* Parameter adjustments */
10337 /* Function Body */
10341 for (i__ = 1; i__ <= i__1; ++i__) {
10342 x += vecte1[i__] * vecte2[i__];
10347 /* ----------------------------------- THE END --------------------------
10353 //=======================================================================
10354 //function : mvcvin2_
10356 //=======================================================================
10357 int mvcvin2_(integer *ncoeff,
10358 doublereal *crvold,
10359 doublereal *crvnew,
10363 /* System generated locals */
10364 integer i__1, i__2;
10366 /* Local variables */
10367 integer m1jm1, ncfm1, j, k;
10369 doublereal cij1, cij2;
10373 /************************************************************************
10378 /* INVERSION OF THE PARAMETERS ON CURVE 2D. */
10382 /* CURVE,2D,INVERSION,PARAMETER. */
10384 /* INPUT ARGUMENTS : */
10385 /* ------------------ */
10386 /* NCOEFF : NB OF COEFF OF THE CURVE. */
10387 /* CRVOLD : CURVE OF ORIGIN */
10389 /* OUTPUT ARGUMENTS : */
10390 /* ------------------- */
10391 /* CRVNEW : THE RESULTING CURVE AFTER CHANGE OF T BY 1-T */
10392 /* IERCOD : 0 OK, */
10393 /* 10 NB OF COEFF NULL OR TOO GREAT. */
10395 /* COMMONS USED : */
10396 /* ---------------- */
10399 /* REFERENCES CALLED : */
10400 /* ---------------------- */
10402 /* DESCRIPTION/NOTES/LIMITATIONS : */
10403 /* ----------------------------------- */
10404 /* THE FOLLOWING CALL IS ABSOLUTELY LEGAL : */
10405 /* CALL MVCVIN2(NCOEFF,CURVE,CURVE,IERCOD), THE TABLE CURVE */
10406 /* BECOMES INPUT AND OUTPUT ARGUMENT (RBD). */
10407 /* BECAUSE OF MCCNP, THE NB OF COEFF OF THE CURVE IS LIMITED TO */
10408 /* NDGCNP+1 = 61. */
10411 /* ***********************************************************************
10415 /* **********************************************************************
10420 /* Serves to provide coefficients of the binome (triangle of Pascal). */
10424 /* Coeff of binome from 0 to 60. read only . init par block data */
10426 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
10427 /* ----------------------------------- */
10428 /* The coefficients of the binome form a triangular matrix. */
10429 /* This matrix is completed in table CNP by transposition. */
10430 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
10432 /* Initialization is done by block-data MMLLL09.RES, */
10433 /* created by program MQINICNP.FOR (see the team (AC) ). */
10437 /* **********************************************************************
10442 /* ***********************************************************************
10445 /* Parameter adjustments */
10449 /* Function Body */
10450 if (*ncoeff < 1 || *ncoeff - 1 > 60) {
10457 /* CONSTANT TERM OF THE NEW CURVE */
10462 for (k = 2; k <= i__1; ++k) {
10463 cij1 += crvold[(k << 1) + 1];
10464 cij2 += crvold[(k << 1) + 2];
10468 if (*ncoeff == 1) {
10472 /* INTERMEDIARY POWERS OF THE PARAMETER */
10474 ncfm1 = *ncoeff - 1;
10477 for (j = 2; j <= i__1; ++j) {
10479 cij1 = crvold[(j << 1) + 1];
10480 cij2 = crvold[(j << 1) + 2];
10482 for (k = j + 1; k <= i__2; ++k) {
10483 bid = mmcmcnp_.cnp[k - 1 + (j - 1) * 61];
10484 cij1 += crvold[(k << 1) + 1] * bid;
10485 cij2 += crvold[(k << 1) + 2] * bid;
10487 crvnew[(j << 1) + 1] = cij1 * m1jm1;
10488 crvnew[(j << 1) + 2] = cij2 * m1jm1;
10491 /* TERM OF THE HIGHEST DEGREE */
10493 crvnew[(*ncoeff << 1) + 1] = -crvold[(*ncoeff << 1) + 1] * m1jm1;
10494 crvnew[(*ncoeff << 1) + 2] = -crvold[(*ncoeff << 1) + 2] * m1jm1;
10498 AdvApp2Var_SysBase::maermsg_("MVCVIN2", iercod, 7L);
10503 //=======================================================================
10504 //function : mvcvinv_
10506 //=======================================================================
10507 int mvcvinv_(integer *ncoeff,
10508 doublereal *crvold,
10509 doublereal *crvnew,
10513 /* System generated locals */
10514 integer i__1, i__2;
10516 /* Local variables */
10517 integer m1jm1, ncfm1, j, k;
10519 //extern /* Subroutine */ int maermsg_();
10520 doublereal cij1, cij2, cij3;
10523 /* **********************************************************************
10528 /* INVERSION OF THE PARAMETER ON A CURBE 3D (I.E. INVERSION */
10529 /* OF THE DIRECTION OF PARSING). */
10533 /* CURVE,INVERSION,PARAMETER. */
10535 /* INPUT ARGUMENTS : */
10536 /* ------------------ */
10537 /* NCOEFF : NB OF COEFF OF THE CURVE. */
10538 /* CRVOLD : CURVE OF ORIGIN */
10540 /* OUTPUT ARGUMENTS : */
10541 /* ------------------- */
10542 /* CRVNEW : RESULTING CURVE AFTER CHANGE OF T INTO 1-T */
10543 /* IERCOD : 0 OK, */
10544 /* 10 NB OF COEFF NULL OR TOO GREAT. */
10546 /* COMMONS USED : */
10547 /* ---------------- */
10550 /* REFERENCES CALLED : */
10551 /* ---------------------- */
10553 /* DESCRIPTION/NOTES/LIMITATIONS : */
10554 /* ----------------------------------- */
10555 /* THE FOLLOWING CALL IS ABSOLUTELY LEGAL : */
10556 /* CALL MVCVINV(NCOEFF,CURVE,CURVE,IERCOD), TABLE CURVE */
10557 /* BECOMES INPUT AND OUTPUT ARGUMENT (RBD). */
10558 /* THE NUMBER OF COEFF OF THE CURVE IS LIMITED TO NDGCNP+1 = 61 */
10559 /* BECAUSE OF USE OF COMMON MCCNP. */
10561 /* ***********************************************************************
10564 /* **********************************************************************
10569 /* Serves to provide the binomial coefficients (triangle of Pascal). */
10573 /* Binomial Coeff from 0 to 60. read only . init par block data */
10575 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
10576 /* ----------------------------------- */
10577 /* The binomial coefficients form a triangular matrix. */
10578 /* This matrix is completed in table CNP by its transposition. */
10579 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
10581 /* Initialisation is done by block-data MMLLL09.RES, */
10582 /* created by program MQINICNP.FOR (see the team (AC) ). */
10584 /* **********************************************************************
10589 /* ***********************************************************************
10592 /* Parameter adjustments */
10596 /* Function Body */
10597 if (*ncoeff < 1 || *ncoeff - 1 > 60) {
10603 /* CONSTANT TERM OF THE NEW CURVE */
10609 for (k = 2; k <= i__1; ++k) {
10610 cij1 += crvold[k * 3 + 1];
10611 cij2 += crvold[k * 3 + 2];
10612 cij3 += crvold[k * 3 + 3];
10618 if (*ncoeff == 1) {
10622 /* INTERMEDIARY POWER OF THE PARAMETER */
10624 ncfm1 = *ncoeff - 1;
10627 for (j = 2; j <= i__1; ++j) {
10629 cij1 = crvold[j * 3 + 1];
10630 cij2 = crvold[j * 3 + 2];
10631 cij3 = crvold[j * 3 + 3];
10633 for (k = j + 1; k <= i__2; ++k) {
10634 bid = mmcmcnp_.cnp[k - 1 + (j - 1) * 61];
10635 cij1 += crvold[k * 3 + 1] * bid;
10636 cij2 += crvold[k * 3 + 2] * bid;
10637 cij3 += crvold[k * 3 + 3] * bid;
10640 crvnew[j * 3 + 1] = cij1 * m1jm1;
10641 crvnew[j * 3 + 2] = cij2 * m1jm1;
10642 crvnew[j * 3 + 3] = cij3 * m1jm1;
10646 /* TERM OF THE HIGHEST DEGREE */
10648 crvnew[*ncoeff * 3 + 1] = -crvold[*ncoeff * 3 + 1] * m1jm1;
10649 crvnew[*ncoeff * 3 + 2] = -crvold[*ncoeff * 3 + 2] * m1jm1;
10650 crvnew[*ncoeff * 3 + 3] = -crvold[*ncoeff * 3 + 3] * m1jm1;
10653 AdvApp2Var_SysBase::maermsg_("MVCVINV", iercod, 7L);
10657 //=======================================================================
10658 //function : mvgaus0_
10660 //=======================================================================
10661 int mvgaus0_(integer *kindic,
10662 doublereal *urootl,
10663 doublereal *hiltab,
10668 /* System generated locals */
10671 /* Local variables */
10672 doublereal tamp[40];
10673 integer ndegl, kg, ii;
10675 /* **********************************************************************
10680 /* Loading of a degree gives roots of LEGENDRE polynom */
10681 /* DEFINED on [-1,1] and weights of Gauss quadrature formulas */
10682 /* (based on corresponding LAGRANGIAN interpolators). */
10683 /* The symmetry relative to 0 is used between [-1,0] and [0,1]. */
10687 /* . VOLUMIC, LEGENDRE, LAGRANGE, GAUSS */
10689 /* INPUT ARGUMENTSE : */
10690 /* ------------------ */
10692 /* KINDIC : Takes values from 1 to 10 depending of the degree */
10693 /* of the used polynom. */
10694 /* The degree of the polynom is equal to 4 k, i.e. 4, 8, */
10695 /* 12, 16, 20, 24, 28, 32, 36 and 40. */
10697 /* OUTPUT ARGUMENTS : */
10698 /* ------------------- */
10700 /* UROOTL : Roots of LEGENDRE polynom in domain [1,0] */
10701 /* given in decreasing order. For domain [-1,0], it is */
10702 /* necessary to take the opposite values. */
10703 /* HILTAB : LAGRANGE interpolators associated to roots. For */
10704 /* opposed roots, interpolatorsare equal. */
10705 /* NBRVAL : Nb of coefficients. Is equal to the half of degree */
10706 /* depending on the symmetry (i.e. 2*KINDIC). */
10708 /* IERCOD : Error code: */
10709 /* < 0 ==> Attention - Warning */
10710 /* =-1 ==> Value of false KINDIC. NBRVAL is forced to 20 */
10712 /* = 0 ==> Everything is OK */
10714 /* COMMON USED : */
10715 /* ---------------- */
10717 /* REFERENCES CALLED : */
10718 /* ------------------- */
10720 /* DESCRIPTION/NOTES/LIMITATIONS : */
10721 /* --------------------------------- */
10722 /* If KINDIC is not correct (i.e < 1 or > 10), the degree is set */
10723 /* to 40 directly (ATTENTION to overload - to avoid it, */
10724 /* preview UROOTL and HILTAB dimensioned at least to 20). */
10726 /* The value of coefficients was calculated with quadruple precision */
10727 /* by JJM with help of GD. */
10728 /* Checking of roots was done by GD. */
10730 /* See detailed explications on the listing */
10732 /* **********************************************************************
10736 /* ------------------------------------ */
10737 /* ****** Test validity of KINDIC ** */
10738 /* ------------------------------------ */
10740 /* Parameter adjustments */
10744 /* Function Body */
10747 if (kg < 1 || kg > 10) {
10752 ndegl = *nbrval << 1;
10754 /* ----------------------------------------------------------------------
10756 /* ****** Load NBRVAL positive roots depending on the degree **
10758 /* ----------------------------------------------------------------------
10760 /* ATTENTION : Sign minus (-) in the loop is intentional. */
10762 mmextrl_(&ndegl, tamp);
10764 for (ii = 1; ii <= i__1; ++ii) {
10765 urootl[ii] = -tamp[ii - 1];
10769 /* ------------------------------------------------------------------- */
10770 /* ****** Loading of NBRVAL Gauss weight depending on the degree ** */
10771 /* ------------------------------------------------------------------- */
10773 mmexthi_(&ndegl, tamp);
10775 for (ii = 1; ii <= i__1; ++ii) {
10776 hiltab[ii] = tamp[ii - 1];
10780 /* ------------------------------- */
10781 /* ****** End of sub-program ** */
10782 /* ------------------------------- */
10787 //=======================================================================
10788 //function : mvpscr2_
10790 //=======================================================================
10791 int mvpscr2_(integer *ncoeff,
10792 doublereal *curve2,
10793 doublereal *tparam,
10794 doublereal *pntcrb)
10796 /* System generated locals */
10799 /* Local variables */
10801 doublereal xxx, yyy;
10805 /* **********************************************************************
10810 /* POSITIONING ON CURVE (NCF,2) IN SPACE OF DIMENSION 2. */
10814 /* TOUS,MATH_ACCES:: COURBE&,POSITIONNEMENT,&POINT. */
10816 /* INPUT ARGUMENTS : */
10817 /* ------------------ */
10818 /* NCOEFF : NUMBER OF COEFFICIENTS OF THE CURVE */
10819 /* CURVE2 : EQUATION OF CURVE 2D */
10820 /* TPARAM : VALUE OF PARAMETER AT GIVEN POINT */
10822 /* OUTPUT ARGUMENTS : */
10823 /* ------------------- */
10824 /* PNTCRB : COORDINATES OF POINT CORRESPONDING TO PARAMETER */
10825 /* TPARAM ON CURVE 2D CURVE2. */
10827 /* COMMONS USED : */
10828 /* ---------------- */
10830 /* REFERENCES CALLED : */
10831 /* ---------------------- */
10833 /* DESCRIPTION/NOTES/LIMITATIONS : */
10834 /* ----------------------------------- */
10835 /* MSCHEMA OF HORNER. */
10838 /* **********************************************************************
10842 /* -------- INITIALIZATIONS AND PROCESSING OF PARTICULAR CASES ----------
10845 /* ---> Cas when NCOEFF > 1 (case STANDARD). */
10846 /* Parameter adjustments */
10850 /* Function Body */
10851 if (*ncoeff >= 2) {
10854 /* ---> Case when NCOEFF <= 1. */
10855 if (*ncoeff <= 0) {
10859 } else if (*ncoeff == 1) {
10860 pntcrb[1] = curve2[3];
10861 pntcrb[2] = curve2[4];
10865 /* -------------------- MSCHEMA OF HORNER (PARTICULAR CASE) --------------
10870 if (*tparam == 1.) {
10874 for (kk = 1; kk <= i__1; ++kk) {
10875 xxx += curve2[(kk << 1) + 1];
10876 yyy += curve2[(kk << 1) + 2];
10880 } else if (*tparam == 0.) {
10881 pntcrb[1] = curve2[3];
10882 pntcrb[2] = curve2[4];
10886 /* ---------------------------- MSCHEMA OF HORNER ------------------------
10888 /* ---> TPARAM is different from 1.D0 and 0.D0. */
10890 ndeg = *ncoeff - 1;
10891 xxx = curve2[(*ncoeff << 1) + 1];
10892 yyy = curve2[(*ncoeff << 1) + 2];
10893 for (kk = ndeg; kk >= 1; --kk) {
10894 xxx = xxx * *tparam + curve2[(kk << 1) + 1];
10895 yyy = yyy * *tparam + curve2[(kk << 1) + 2];
10900 /* ------------------------ RECOVER THE CALCULATED POINT ---------------
10907 /* ------------------------------ THE END -------------------------------
10914 //=======================================================================
10915 //function : mvpscr3_
10917 //=======================================================================
10918 int mvpscr3_(integer *ncoeff,
10919 doublereal *curve3,
10920 doublereal *tparam,
10921 doublereal *pntcrb)
10924 /* System generated locals */
10927 /* Local variables */
10929 doublereal xxx, yyy, zzz;
10933 /* **********************************************************************
10938 /* POSITIONING ON A CURVE (3,NCF) IN THE SPACE OF DIMENSION 3. */
10942 /* TOUS, MATH_ACCES:: COURBE&,POSITIONNEMENT,&POINT. */
10944 /* INPUT ARGUMENTS : */
10945 /* ------------------ */
10946 /* NCOEFF : NB OF COEFFICIENTS OF THE CURVE */
10947 /* CURVE3 : EQUATION OF CURVE 3D */
10948 /* TPARAM : VALUE OF THE PARAMETER AT THE GIVEN POINT */
10950 /* OUTPUT ARGUMENTS : */
10951 /* ------------------- */
10952 /* PNTCRB : COORDINATES OF THE POINT CORRESPONDING TO PARAMETER */
10953 /* TPARAM ON CURVE 3D CURVE3. */
10955 /* COMMONS USED : */
10956 /* ---------------- */
10958 /* REFERENCES CALLED : */
10959 /* ---------------------- */
10962 /* DESCRIPTION/NOTES/LIMITATIONS : */
10963 /* ----------------------------------- */
10964 /* MSCHEMA OF HORNER. */
10966 /* **********************************************************************
10969 /* **********************************************************************
10973 /* -------- INITIALISATIONS AND PROCESSING OF PARTICULAR CASES ----------
10976 /* ---> Case when NCOEFF > 1 (cas STANDARD). */
10977 /* Parameter adjustments */
10981 /* Function Body */
10982 if (*ncoeff >= 2) {
10985 /* ---> Case when NCOEFF <= 1. */
10986 if (*ncoeff <= 0) {
10991 } else if (*ncoeff == 1) {
10992 pntcrb[1] = curve3[4];
10993 pntcrb[2] = curve3[5];
10994 pntcrb[3] = curve3[6];
10998 /* -------------------- MSCHEMA OF HORNER (PARTICULAR CASE) --------------
11003 if (*tparam == 1.) {
11008 for (kk = 1; kk <= i__1; ++kk) {
11009 xxx += curve3[kk * 3 + 1];
11010 yyy += curve3[kk * 3 + 2];
11011 zzz += curve3[kk * 3 + 3];
11015 } else if (*tparam == 0.) {
11016 pntcrb[1] = curve3[4];
11017 pntcrb[2] = curve3[5];
11018 pntcrb[3] = curve3[6];
11022 /* ---------------------------- MSCHEMA OF HORNER ------------------------
11024 /* ---> Here TPARAM is different from 1.D0 and 0.D0. */
11026 ndeg = *ncoeff - 1;
11027 xxx = curve3[*ncoeff * 3 + 1];
11028 yyy = curve3[*ncoeff * 3 + 2];
11029 zzz = curve3[*ncoeff * 3 + 3];
11030 for (kk = ndeg; kk >= 1; --kk) {
11031 xxx = xxx * *tparam + curve3[kk * 3 + 1];
11032 yyy = yyy * *tparam + curve3[kk * 3 + 2];
11033 zzz = zzz * *tparam + curve3[kk * 3 + 3];
11038 /* ------------------------ RETURN THE CALCULATED POINT ------------------
11046 /* ------------------------------ THE END -------------------------------
11053 //=======================================================================
11054 //function : AdvApp2Var_MathBase::mvsheld_
11056 //=======================================================================
11057 int AdvApp2Var_MathBase::mvsheld_(integer *n,
11063 /* System generated locals */
11064 integer dtab_dim1, dtab_offset, i__1, i__2;
11066 /* Local variables */
11069 integer i3, i4, i5, incrp1;
11072 /************************************************************************
11077 /* PARSING OF COLUMNS OF TABLE OF REAL*8 BY SHELL METHOD*/
11078 /* (IN INCREASING ORDER) */
11082 /* POINT-ENTRY, PARSING, SHELL */
11084 /* INPUT ARGUMENTS : */
11085 /* ------------------ */
11086 /* N : NUMBER OF COLUMNS OF THE TABLE */
11087 /* IS : NUMBER OF LINE OF THE TABLE */
11088 /* DTAB : TABLE OF REAL*8 TO BE PARSED */
11089 /* ICLE : POSITION OF THE KEY ON THE COLUMN */
11091 /* OUTPUT ARGUMENTS : */
11092 /* ------------------- */
11093 /* DTAB : PARSED TABLE */
11095 /* COMMONS USED : */
11096 /* ---------------- */
11099 /* REFERENCES CALLED : */
11100 /* ---------------------- */
11103 /* DESCRIPTION/NOTES/LIMITATIONS : */
11104 /* ----------------------------------- */
11105 /* CLASSIC SHELL METHOD : PARSING BY SERIES */
11106 /* Declaration DTAB(IS, 1) corresponds to DTAB(IS, *) */
11108 /* ***********************************************************************
11112 /* Parameter adjustments */
11114 dtab_offset = dtab_dim1 + 1;
11115 dtab -= dtab_offset;
11117 /* Function Body */
11121 /* ------------------------ */
11123 /* INITIALIZATION OF THE SEQUENCE OF INCREMENTS */
11124 /* FIND THE GREATEST INCREMENT SO THAT INCR < N/9 */
11128 if (incr >= *n / 9) {
11131 /* ----------------------------- */
11132 incr = incr * 3 + 1;
11135 /* LOOP ON INCREMENTS TILL INCR = 1 */
11136 /* PARSING BY SERIES DISTANT FROM INCR */
11140 /* ----------------- */
11142 for (i3 = incrp1; i3 <= i__1; ++i3) {
11143 /* ---------------------- */
11145 /* SET ELEMENT I3 AT ITS PLACE IN THE SERIES */
11152 /* ------------------------- */
11153 if (dtab[*icle + i4 * dtab_dim1] <= dtab[*icle + (i4 + incr) *
11159 for (i5 = 1; i5 <= i__2; ++i5) {
11160 /* ------------------ */
11161 dsave = dtab[i5 + i4 * dtab_dim1];
11162 dtab[i5 + i4 * dtab_dim1] = dtab[i5 + (i4 + incr) * dtab_dim1];
11163 dtab[i5 + (i4 + incr) * dtab_dim1] = dsave;
11174 /* PASSAGE TO THE NEXT INCREMENT */
11185 //=======================================================================
11186 //function : AdvApp2Var_MathBase::mzsnorm_
11188 //=======================================================================
11189 doublereal AdvApp2Var_MathBase::mzsnorm_(integer *ndimen,
11190 doublereal *vecteu)
11193 /* System generated locals */
11195 doublereal ret_val, d__1, d__2;
11197 /* Local variables */
11199 integer i__, irmax;
11203 /* ***********************************************************************
11208 /* SERVES to calculate the euclidian norm of a vector : */
11209 /* ____________________________ */
11210 /* Z = V V(1)**2 + V(2)**2 + ... */
11216 /* INPUT ARGUMENTS : */
11217 /* ------------------ */
11218 /* NDIMEN : Dimension of the vector */
11219 /* VECTEU : vector of dimension NDIMEN */
11221 /* OUTPUT ARGUMENTS : */
11222 /* ------------------- */
11223 /* MZSNORM : Value of the euclidian norm of vector VECTEU */
11225 /* COMMONS USED : */
11226 /* ---------------- */
11230 /* REFERENCES CALLED : */
11231 /* ---------------------- */
11233 /* R*8 ABS R*8 SQRT */
11235 /* DESCRIPTION/NOTESS/LIMITATIONS : */
11236 /* ----------------------------------- */
11237 /* To limit the risks of overflow, */
11238 /* the term of the strongest absolute value is factorized : */
11239 /* _______________________ */
11240 /* Z = !V(1)! * V 1 + (V(2)/V(1))**2 + ... */
11243 /* ***********************************************************************
11246 /* ***********************************************************************
11250 /* ***********************************************************************
11253 /* ***********************************************************************
11256 /* ___ Find the strongest absolute value term */
11258 /* Parameter adjustments */
11261 /* Function Body */
11264 for (i__ = 2; i__ <= i__1; ++i__) {
11265 if ((d__1 = vecteu[irmax], advapp_abs(d__1)) < (d__2 = vecteu[i__], advapp_abs(d__2)
11272 /* ___ Calculate the norme */
11274 if ((d__1 = vecteu[irmax], advapp_abs(d__1)) < 1.) {
11277 for (i__ = 1; i__ <= i__1; ++i__) {
11278 /* Computing 2nd power */
11279 d__1 = vecteu[i__];
11280 xsom += d__1 * d__1;
11283 ret_val = sqrt(xsom);
11287 for (i__ = 1; i__ <= i__1; ++i__) {
11288 if (i__ == irmax) {
11291 /* Computing 2nd power */
11292 d__1 = vecteu[i__] / vecteu[irmax];
11293 xsom += d__1 * d__1;
11297 ret_val = (d__1 = vecteu[irmax], advapp_abs(d__1)) * sqrt(xsom);
11300 /* ***********************************************************************
11302 /* RETURN CALLING PROGRAM */
11303 /* ***********************************************************************