2 // AdvApp2Var_MathBase.cxx
5 #include <AdvApp2Var_SysBase.hxx>
6 #include <AdvApp2Var_Data_f2c.hxx>
7 #include <AdvApp2Var_MathBase.hxx>
8 #include <AdvApp2Var_Data.hxx>
12 int mmchole_(integer *mxcoef,
24 int mmrslss_(integer *mxcoef,
34 int mfac_(doublereal *f,
38 int mmaper0_(integer *ncofmx,
46 int mmaper2_(integer *ncofmx,
55 int mmaper4_(integer *ncofmx,
64 int mmaper6_(integer *ncofmx,
73 int mmarc41_(integer *ndimax,
83 int mmatvec_(integer *nligne,
94 int mmcvstd_(integer *ncofmx,
102 int mmdrvcb_(integer *ideriv,
111 int mmexthi_(integer *ndegre,
115 int mmextrl_(integer *ndegre,
121 int mmherm0_(doublereal *debfin,
125 int mmherm1_(doublereal *debfin,
131 int mmloncv_(integer *ndimax,
140 int mmpojac_(doublereal *tparam,
148 int mmrslw_(integer *normax,
156 int mmtmave_(integer *nligne,
165 int mmtrpj0_(integer *ncofmx,
174 int mmtrpj2_(integer *ncofmx,
184 int mmtrpj4_(integer *ncofmx,
193 int mmtrpj6_(integer *ncofmx,
202 integer pow__ii(integer *x,
206 int mvcvin2_(integer *ncoeff,
212 int mvcvinv_(integer *ncoeff,
218 int mvgaus0_(integer *kindic,
224 int mvpscr2_(integer *ncoeff,
230 int mvpscr3_(integer *ncoeff,
236 doublereal eps1, eps2, eps3, eps4;
237 integer niterm, niterr;
241 doublereal tdebut, tfinal, verifi, cmherm[576];
244 //=======================================================================
245 //function : AdvApp2Var_MathBase::mdsptpt_
247 //=======================================================================
248 int AdvApp2Var_MathBase::mdsptpt_(integer *ndimen,
254 static integer c__8 = 8;
255 /* System generated locals */
259 /* Local variables */
261 static doublereal differ[100];
265 /* **********************************************************************
270 /* CALCULATE DISTANCE BETWEEN TWO POINTS */
274 /* DISTANCE,POINT. */
276 /* INPUT ARGUMENTS : */
277 /* ------------------ */
278 /* NDIMEN: Space Dimension. */
279 /* POINT1: Table of coordinates of the 1st point. */
280 /* POINT2: Table of coordinates of the 2nd point. */
282 /* OUTPUT ARGUMENTS : */
283 /* ------------------- */
284 /* DISTAN: Distance between 2 points. */
287 /* ---------------- */
289 /* REFERENCES CALLED : */
290 /* ----------------------- */
292 /* DESCRIPTION/NOTES/LIMITATIONS : */
293 /* ----------------------------------- */
295 /* **********************************************************************
299 /* ***********************************************************************
302 /* ***********************************************************************
305 /* Parameter adjustment */
313 /* ***********************************************************************
316 /* ***********************************************************************
320 AdvApp2Var_SysBase::mcrrqst_(&c__8, ndimen, differ, &iofset, &ier);
323 /* --- If allocation is refused, the trivial method is applied. */
329 for (i__ = 1; i__ <= i__1; ++i__) {
330 /* Computing 2nd power */
331 d__1 = point1[i__] - point2[i__];
332 *distan += d__1 * d__1;
334 *distan = sqrt(*distan);
336 /* --- Otherwise MZSNORM is used to minimize the risks of overflow
341 for (i__ = 1; i__ <= i__1; ++i__) {
343 differ[j] = point2[i__] - point1[i__];
346 *distan = AdvApp2Var_MathBase::mzsnorm_(ndimen, &differ[iofset]);
350 /* ***********************************************************************
352 /* RETURN CALLING PROGRAM */
353 /* ***********************************************************************
356 /* --- Dynamic Desallocation */
359 AdvApp2Var_SysBase::mcrdelt_(&c__8, ndimen, differ, &iofset, &ier);
365 //=======================================================================
368 //=======================================================================
369 int mfac_(doublereal *f,
373 /* System generated locals */
376 /* Local variables */
379 /* FORTRAN CONFORME AU TEXT */
380 /* CALCUL DE MFACTORIEL N */
381 /* Parameter adjustments */
387 for (i__ = 2; i__ <= i__1; ++i__) {
389 f[i__] = i__ * f[i__ - 1];
394 //=======================================================================
395 //function : AdvApp2Var_MathBase::mmapcmp_
397 //=======================================================================
398 int AdvApp2Var_MathBase::mmapcmp_(integer *ndim,
405 /* System generated locals */
406 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
409 /* Local variables */
410 static integer ipair, nd, ndegre, impair, ibb, idg;
411 //extern int mgsomsg_();//mgenmsg_(),
415 /* **********************************************************************
420 /* Compression of curve CRVOLD in a table of */
421 /* coeff. of even : CRVNEW(*,0,*) */
422 /* and uneven range : CRVNEW(*,1,*). */
426 /* COMPRESSION,CURVE. */
428 /* INPUT ARGUMENTS : */
429 /* ------------------ */
430 /* NDIM : Space Dimension. */
431 /* NCOFMX : Max nb of coeff. of the curve to compress. */
432 /* NCOEFF : Max nb of coeff. of the compressed curve. */
433 /* CRVOLD : The curve (0:NCOFMX-1,NDIM) to compress. */
435 /* OUTPUT ARGUMENTS : */
436 /* ------------------- */
437 /* CRVNEW : Curve compacted in (0:(NCOEFF-1)/2,0,NDIM) (containing
439 /* even terms) and in (0:(NCOEFF-1)/2,1,NDIM) */
440 /* (containing uneven terms). */
443 /* ---------------- */
445 /* REFERENCES CALLED : */
446 /* ----------------------- */
448 /* DESCRIPTION/NOTES/LIMITATIONS : */
449 /* ----------------------------------- */
450 /* This routine is useful to prepare coefficients of a */
451 /* curve in an orthogonal base (Legendre or Jacobi) before */
452 /* calculating the coefficients in the canonical; base [-1,1] by */
454 /* ***********************************************************************
457 /* Name of the routine */
459 /* Parameter adjustments */
460 crvold_dim1 = *ncofmx;
461 crvold_offset = crvold_dim1;
462 crvold -= crvold_offset;
463 crvnew_dim1 = (*ncoeff - 1) / 2 + 1;
464 crvnew_offset = crvnew_dim1 << 1;
465 crvnew -= crvnew_offset;
468 ibb = AdvApp2Var_SysBase::mnfndeb_();
470 AdvApp2Var_SysBase::mgenmsg_("MMAPCMP", 7L);
473 ndegre = *ncoeff - 1;
475 for (nd = 1; nd <= i__1; ++nd) {
478 for (idg = 0; idg <= i__2; ++idg) {
479 crvnew[idg + (nd << 1) * crvnew_dim1] = crvold[ipair + nd *
488 i__2 = (ndegre - 1) / 2;
489 for (idg = 0; idg <= i__2; ++idg) {
490 crvnew[idg + ((nd << 1) + 1) * crvnew_dim1] = crvold[impair + nd *
501 /* ---------------------------------- The end ---------------------------
505 AdvApp2Var_SysBase::mgsomsg_("MMAPCMP", 7L);
510 //=======================================================================
511 //function : mmaper0_
513 //=======================================================================
514 int mmaper0_(integer *ncofmx,
523 /* System generated locals */
524 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
527 /* Local variables */
529 static doublereal bidon;
530 static integer ii, nd;
533 /* ***********************************************************************
538 /* Calculate the max error of approximation done when */
539 /* only the first NCFNEW coefficients of a curve are preserved.
541 /* Degree NCOEFF-1 written in the base of Legendre (Jacobi */
546 /* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
548 /* INPUT ARGUMENTS : */
549 /* ------------------ */
550 /* NCOFMX : Max. degree of the curve. */
551 /* NDIMEN : Space dimension. */
552 /* NCOEFF : Degree +1 of the curve. */
553 /* CRVLGD : Curve the degree which of should be lowered. */
554 /* NCFNEW : Degree +1 of the resulting polynom. */
556 /* OUTPUT ARGUMENTS : */
557 /* ------------------- */
558 /* YCVMAX : Auxiliary Table (max error on each dimension).
560 /* ERRMAX : Precision of the approximation. */
563 /* ---------------- */
565 /* REFERENCES CALLED : */
566 /* ----------------------- */
568 /* DESCRIPTION/NOTES/LIMITATIONS : */
569 /* ----------------------------------- */
570 /* ***********************************************************************
574 /* ------------------- Init to calculate an error -----------------------
577 /* Parameter adjustments */
579 crvlgd_dim1 = *ncofmx;
580 crvlgd_offset = crvlgd_dim1 + 1;
581 crvlgd -= crvlgd_offset;
585 for (ii = 1; ii <= i__1; ++ii) {
590 /* ------ Minimum that can be reached : Stop at 1 or NCFNEW ------
594 if (*ncfnew + 1 > ncut) {
598 /* -------------- Elimination of high degree coefficients-----------
600 /* ----------- Loop on the series of Legendre: NCUT --> NCOEFF --------
604 for (ii = ncut; ii <= i__1; ++ii) {
605 /* Factor of renormalization (Maximum of Li(t)). */
606 bidon = ((ii - 1) * 2. + 1.) / 2.;
610 for (nd = 1; nd <= i__2; ++nd) {
611 ycvmax[nd] += (d__1 = crvlgd[ii + nd * crvlgd_dim1], abs(d__1)) *
618 /* -------------- The error is the norm of the vector error ---------------
621 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
623 /* --------------------------------- Fin --------------------------------
629 //=======================================================================
630 //function : mmaper2_
632 //=======================================================================
633 int mmaper2_(integer *ncofmx,
642 /* Initialized data */
644 static doublereal xmaxj[57] = { .9682458365518542212948163499456,
645 .986013297183269340427888048593603,
646 1.07810420343739860362585159028115,
647 1.17325804490920057010925920756025,
648 1.26476561266905634732910520370741,
649 1.35169950227289626684434056681946,
650 1.43424378958284137759129885012494,
651 1.51281316274895465689402798226634,
652 1.5878364329591908800533936587012,
653 1.65970112228228167018443636171226,
654 1.72874345388622461848433443013543,
655 1.7952515611463877544077632304216,
656 1.85947199025328260370244491818047,
657 1.92161634324190018916351663207101,
658 1.98186713586472025397859895825157,
659 2.04038269834980146276967984252188,
660 2.09730119173852573441223706382076,
661 2.15274387655763462685970799663412,
662 2.20681777186342079455059961912859,
663 2.25961782459354604684402726624239,
664 2.31122868752403808176824020121524,
665 2.36172618435386566570998793688131,
666 2.41117852396114589446497298177554,
667 2.45964731268663657873849811095449,
668 2.50718840313973523778244737914028,
669 2.55385260994795361951813645784034,
670 2.59968631659221867834697883938297,
671 2.64473199258285846332860663371298,
672 2.68902863641518586789566216064557,
673 2.73261215675199397407027673053895,
674 2.77551570192374483822124304745691,
675 2.8177699459714315371037628127545,
676 2.85940333797200948896046563785957,
677 2.90044232019793636101516293333324,
678 2.94091151970640874812265419871976,
679 2.98083391718088702956696303389061,
680 3.02023099621926980436221568258656,
681 3.05912287574998661724731962377847,
682 3.09752842783622025614245706196447,
683 3.13546538278134559341444834866301,
684 3.17295042316122606504398054547289,
685 3.2099992681699613513775259670214,
686 3.24662674946606137764916854570219,
687 3.28284687953866689817670991319787,
688 3.31867291347259485044591136879087,
689 3.35411740487202127264475726990106,
690 3.38919225660177218727305224515862,
691 3.42390876691942143189170489271753,
692 3.45827767149820230182596660024454,
693 3.49230918177808483937957161007792,
694 3.5260130200285724149540352829756,
695 3.55939845146044235497103883695448,
696 3.59247431368364585025958062194665,
697 3.62524904377393592090180712976368,
698 3.65773070318071087226169680450936,
699 3.68992700068237648299565823810245,
700 3.72184531357268220291630708234186 };
702 /* System generated locals */
703 integer crvjac_dim1, crvjac_offset, i__1, i__2;
706 /* Local variables */
707 static integer idec, ncut;
708 static doublereal bidon;
709 static integer ii, nd;
713 /* ***********************************************************************
718 /* Calculate max approximation error i faite lorsque l' on */
719 /* ne conserve que les premiers NCFNEW coefficients d' une courbe
721 /* de degre NCOEFF-1 ecrite dans la base de Jacobi d' ordre 2. */
725 /* JACOBI, POLYGON, APPROXIMATION, ERROR. */
727 /* INPUT ARGUMENTS : */
728 /* ------------------ */
729 /* NCOFMX : Max. degree of the curve. */
730 /* NDIMEN : Space dimension. */
731 /* NCOEFF : Degree +1 of the curve. */
732 /* CRVLGD : Curve the degree which of should be lowered. */
733 /* NCFNEW : Degree +1 of the resulting polynom. */
735 /* OUTPUT ARGUMENTS : */
736 /* ------------------- */
737 /* YCVMAX : Auxiliary Table (max error on each dimension).
739 /* ERRMAX : Precision of the approximation. */
742 /* ---------------- */
744 /* REFERENCES CALLED : */
745 /* ----------------------- */
746 /* DESCRIPTION/NOTES/LIMITATIONS : */
747 /* ----------------------------------- */
751 /* ------------------ Table of maximums of (1-t2)*Ji(t) ----------------
754 /* Parameter adjustments */
756 crvjac_dim1 = *ncofmx;
757 crvjac_offset = crvjac_dim1 + 1;
758 crvjac -= crvjac_offset;
764 /* ------------------- Init for error calculation -----------------------
768 for (ii = 1; ii <= i__1; ++ii) {
773 /* ------ Min. Degree that can be attained : Stop at 3 or NCFNEW ------
778 i__1 = idec, i__2 = *ncfnew + 1;
779 ncut = max(i__1,i__2);
781 /* -------------- Removal of coefficients of high degree -----------
783 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
787 for (ii = ncut; ii <= i__1; ++ii) {
788 /* Factor of renormalization. */
789 bidon = xmaxj[ii - idec];
791 for (nd = 1; nd <= i__2; ++nd) {
792 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], abs(d__1)) *
799 /* -------------- The error is the norm of the vector error ---------------
802 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
804 /* --------------------------------- Fin --------------------------------
810 /* MAPER4.f -- translated by f2c (version 19960827).
811 You must link the resulting object file with the libraries:
812 -lf2c -lm (in that order)
816 //=======================================================================
817 //function : mmaper4_
819 //=======================================================================
820 int mmaper4_(integer *ncofmx,
828 /* Initialized data */
830 static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
831 1.05299572648705464724876659688996,
832 1.0949715351434178709281698645813,
833 1.15078388379719068145021100764647,
834 1.2094863084718701596278219811869,
835 1.26806623151369531323304177532868,
836 1.32549784426476978866302826176202,
837 1.38142537365039019558329304432581,
838 1.43575531950773585146867625840552,
839 1.48850442653629641402403231015299,
840 1.53973611681876234549146350844736,
841 1.58953193485272191557448229046492,
842 1.63797820416306624705258190017418,
843 1.68515974143594899185621942934906,
844 1.73115699602477936547107755854868,
845 1.77604489805513552087086912113251,
846 1.81989256661534438347398400420601,
847 1.86276344480103110090865609776681,
848 1.90471563564740808542244678597105,
849 1.94580231994751044968731427898046,
850 1.98607219357764450634552790950067,
851 2.02556989246317857340333585562678,
852 2.06433638992049685189059517340452,
853 2.10240936014742726236706004607473,
854 2.13982350649113222745523925190532,
855 2.17661085564771614285379929798896,
856 2.21280102016879766322589373557048,
857 2.2484214321456956597803794333791,
858 2.28349755104077956674135810027654,
859 2.31805304852593774867640120860446,
860 2.35210997297725685169643559615022,
861 2.38568889602346315560143377261814,
862 2.41880904328694215730192284109322,
863 2.45148841120796359750021227795539,
864 2.48374387161372199992570528025315,
865 2.5155912654873773953959098501893,
866 2.54704548720896557684101746505398,
867 2.57812056037881628390134077704127,
868 2.60882970619319538196517982945269,
869 2.63918540521920497868347679257107,
870 2.66919945330942891495458446613851,
871 2.69888301230439621709803756505788,
872 2.72824665609081486737132853370048,
873 2.75730041251405791603760003778285,
874 2.78605380158311346185098508516203,
875 2.81451587035387403267676338931454,
876 2.84269522483114290814009184272637,
877 2.87060005919012917988363332454033,
878 2.89823818258367657739520912946934,
879 2.92561704377132528239806135133273,
880 2.95274375377994262301217318010209,
881 2.97962510678256471794289060402033,
882 3.00626759936182712291041810228171,
883 3.03267744830655121818899164295959,
884 3.05886060707437081434964933864149 };
886 /* System generated locals */
887 integer crvjac_dim1, crvjac_offset, i__1, i__2;
890 /* Local variables */
891 static integer idec, ncut;
892 static doublereal bidon;
893 static integer ii, nd;
897 /* ***********************************************************************
902 /* Calculate the max. error of approximation made when */
903 /* only first NCFNEW coefficients of a curve are preserved
905 /* degree NCOEFF-1 is written in the base of Jacobi of order 4. */
908 /* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
910 /* INPUT ARGUMENTS : */
911 /* ------------------ */
912 /* NCOFMX : Max. degree of the curve. */
913 /* NDIMEN : Space dimension. */
914 /* NCOEFF : Degree +1 of the curve. */
915 /* CRVJAC : Curve the degree which of should be lowered. */
916 /* NCFNEW : Degree +1 of the resulting polynom. */
918 /* OUTPUT ARGUMENTS : */
919 /* ------------------- */
920 /* YCVMAX : Auxiliary Table (max error on each dimension).
922 /* ERRMAX : Precision of the approximation. */
925 /* ---------------- */
927 /* REFERENCES CALLED : */
928 /* ----------------------- */
930 /* DESCRIPTION/NOTES/LIMITATIONS : */
933 /* ***********************************************************************
937 /* ---------------- Table of maximums of ((1-t2)2)*Ji(t) ---------------
940 /* Parameter adjustments */
942 crvjac_dim1 = *ncofmx;
943 crvjac_offset = crvjac_dim1 + 1;
944 crvjac -= crvjac_offset;
950 /* ------------------- Init for error calculation -----------------------
954 for (ii = 1; ii <= i__1; ++ii) {
959 /* ------ Min. Degree that can be attained : Stop at 5 or NCFNEW ------
964 i__1 = idec, i__2 = *ncfnew + 1;
965 ncut = max(i__1,i__2);
967 /* -------------- Removal of high degree coefficients -----------
969 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
973 for (ii = ncut; ii <= i__1; ++ii) {
974 /* Factor of renormalisation. */
975 bidon = xmaxj[ii - idec];
977 for (nd = 1; nd <= i__2; ++nd) {
978 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], abs(d__1)) *
985 /* -------------- The error is the norm of the error vector ---------------
988 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
990 /* --------------------------------- End --------------------------------
996 //=======================================================================
997 //function : mmaper6_
999 //=======================================================================
1000 int mmaper6_(integer *ncofmx,
1009 /* Initialized data */
1011 static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
1012 1.11626917091567929907256116528817,
1013 1.1327140810290884106278510474203,
1014 1.1679452722668028753522098022171,
1015 1.20910611986279066645602153641334,
1016 1.25228283758701572089625983127043,
1017 1.29591971597287895911380446311508,
1018 1.3393138157481884258308028584917,
1019 1.3821288728999671920677617491385,
1020 1.42420414683357356104823573391816,
1021 1.46546895108549501306970087318319,
1022 1.50590085198398789708599726315869,
1023 1.54550385142820987194251585145013,
1024 1.58429644271680300005206185490937,
1025 1.62230484071440103826322971668038,
1026 1.65955905239130512405565733793667,
1027 1.69609056468292429853775667485212,
1028 1.73193098017228915881592458573809,
1029 1.7671112206990325429863426635397,
1030 1.80166107681586964987277458875667,
1031 1.83560897003644959204940535551721,
1032 1.86898184653271388435058371983316,
1033 1.90180515174518670797686768515502,
1034 1.93410285411785808749237200054739,
1035 1.96589749778987993293150856865539,
1036 1.99721027139062501070081653790635,
1037 2.02806108474738744005306947877164,
1038 2.05846864831762572089033752595401,
1039 2.08845055210580131460156962214748,
1040 2.11802334209486194329576724042253,
1041 2.14720259305166593214642386780469,
1042 2.17600297710595096918495785742803,
1043 2.20443832785205516555772788192013,
1044 2.2325216999457379530416998244706,
1045 2.2602654243075083168599953074345,
1046 2.28768115912702794202525264301585,
1047 2.3147799369092684021274946755348,
1048 2.34157220782483457076721300512406,
1049 2.36806787963276257263034969490066,
1050 2.39427635443992520016789041085844,
1051 2.42020656255081863955040620243062,
1052 2.44586699364757383088888037359254,
1053 2.47126572552427660024678584642791,
1054 2.49641045058324178349347438430311,
1055 2.52130850028451113942299097584818,
1056 2.54596686772399937214920135190177,
1057 2.5703922285006754089328998222275,
1058 2.59459096001908861492582631591134,
1059 2.61856915936049852435394597597773,
1060 2.64233265984385295286445444361827,
1061 2.66588704638685848486056711408168,
1062 2.68923766976735295746679957665724,
1063 2.71238965987606292679677228666411 };
1065 /* System generated locals */
1066 integer crvjac_dim1, crvjac_offset, i__1, i__2;
1069 /* Local variables */
1070 static integer idec, ncut;
1071 static doublereal bidon;
1072 static integer ii, nd;
1076 /* ***********************************************************************
1080 /* Calculate the max. error of approximation made when */
1081 /* only first NCFNEW coefficients of a curve are preserved
1083 /* degree NCOEFF-1 is written in the base of Jacobi of order 6. */
1086 /* JACOBI,POLYGON,APPROXIMATION,ERROR. */
1088 /* INPUT ARGUMENTS : */
1089 /* ------------------ */
1090 /* NCOFMX : Max. degree of the curve. */
1091 /* NDIMEN : Space dimension. */
1092 /* NCOEFF : Degree +1 of the curve. */
1093 /* CRVJAC : Curve the degree which of should be lowered. */
1094 /* NCFNEW : Degree +1 of the resulting polynom. */
1096 /* OUTPUT ARGUMENTS : */
1097 /* ------------------- */
1098 /* YCVMAX : Auxiliary Table (max error on each dimension).
1100 /* ERRMAX : Precision of the approximation. */
1102 /* COMMONS USED : */
1103 /* ---------------- */
1105 /* REFERENCES CALLED : */
1106 /* ----------------------- */
1108 /* DESCRIPTION/NOTES/LIMITATIONS : */
1110 /* ***********************************************************************
1114 /* ---------------- Table of maximums of ((1-t2)3)*Ji(t) ---------------
1117 /* Parameter adjustments */
1119 crvjac_dim1 = *ncofmx;
1120 crvjac_offset = crvjac_dim1 + 1;
1121 crvjac -= crvjac_offset;
1127 /* ------------------- Init for error calculation -----------------------
1131 for (ii = 1; ii <= i__1; ++ii) {
1136 /* ------ Min Degree that can be attained : Stop at 3 or NCFNEW ------
1141 i__1 = idec, i__2 = *ncfnew + 1;
1142 ncut = max(i__1,i__2);
1144 /* -------------- Removal of high degree coefficients -----------
1146 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
1150 for (ii = ncut; ii <= i__1; ++ii) {
1151 /* Factor of renormalization. */
1152 bidon = xmaxj[ii - idec];
1154 for (nd = 1; nd <= i__2; ++nd) {
1155 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], abs(d__1)) *
1162 /* -------------- The error is the norm of the vector error ---------------
1165 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
1167 /* --------------------------------- END --------------------------------
1173 //=======================================================================
1174 //function : AdvApp2Var_MathBase::mmaperx_
1176 //=======================================================================
1177 int AdvApp2Var_MathBase::mmaperx_(integer *ncofmx,
1188 /* System generated locals */
1189 integer crvjac_dim1, crvjac_offset;
1191 /* Local variables */
1192 static integer jord;
1195 /* **********************************************************************
1199 /* Calculate the max. error of approximation made when */
1200 /* only first NCFNEW coefficients of a curve are preserved
1202 /* degree NCOEFF-1 is written in the base of Jacobi of order IORDRE. */
1205 /* JACOBI,LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
1207 /* INPUT ARGUMENTS : */
1208 /* ------------------ */
1209 /* NCOFMX : Max. degree of the curve. */
1210 /* NDIMEN : Space dimension. */
1211 /* NCOEFF : Degree +1 of the curve. */
1212 /* IORDRE : Order of continuity at the extremities. */
1213 /* CRVJAC : Curve the degree which of should be lowered. */
1214 /* NCFNEW : Degree +1 of the resulting polynom. */
1216 /* OUTPUT ARGUMENTS : */
1217 /* ------------------- */
1218 /* YCVMAX : Auxiliary Table (max error on each dimension).
1220 /* ERRMAX : Precision of the approximation. */
1221 /* IERCOD = 0, OK */
1222 /* = 1, order of constraints (IORDRE) is not within the */
1223 /* autorized values. */
1224 /* COMMONS USED : */
1225 /* ---------------- */
1227 /* REFERENCES CALLED : */
1228 /* ----------------------- */
1230 /* DESCRIPTION/NOTES/LIMITATIONS : */
1231 /* ----------------------------------- */
1232 /* Canceled and replaced MMAPERR. */
1233 /* ***********************************************************************
1237 /* Parameter adjustments */
1239 crvjac_dim1 = *ncofmx;
1240 crvjac_offset = crvjac_dim1 + 1;
1241 crvjac -= crvjac_offset;
1245 /* --> Order of Jacobi polynoms */
1246 jord = ( *iordre + 1) << 1;
1249 mmaper0_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1251 } else if (jord == 2) {
1252 mmaper2_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1254 } else if (jord == 4) {
1255 mmaper4_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1257 } else if (jord == 6) {
1258 mmaper6_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1264 /* ----------------------------------- Fin ------------------------------
1270 //=======================================================================
1271 //function : mmarc41_
1273 //=======================================================================
1274 int mmarc41_(integer *ndimax,
1284 /* System generated locals */
1285 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
1288 /* Local variables */
1289 static integer nboct;
1290 static doublereal tbaux[61];
1292 static doublereal bid;
1293 static integer ncf, ncj;
1296 /* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
1297 /* IMPLICIT INTEGER (I-N) */
1299 /* ***********************************************************************
1304 /* Creation of curve C2(v) defined on (0,1) identic to */
1305 /* curve C1(u) defined on (U0,U1) (change of parameter */
1310 /* LIMITATION, RESTRICTION, CURVE */
1312 /* INPUT ARGUMENTS : */
1313 /* ------------------ */
1314 /* NDIMAX : Space Dimensioning. */
1315 /* NDIMEN : Curve Dimension. */
1316 /* NCOEFF : Nb of coefficients of the curve. */
1317 /* CRVOLD : Curve to be limited. */
1318 /* UPARA0 : Min limit of the interval limiting the curve.
1320 /* UPARA1 : Max limit of the interval limiting the curve.
1323 /* OUTPUT ARGUMENTS : */
1324 /* ------------------- */
1325 /* CRVNEW : Relimited curve, defined on (0,1) and equal to */
1326 /* CRVOLD defined on (U0,U1). */
1327 /* IERCOD : = 0, OK */
1328 /* =10, Nb of coeff. <1 or > 61. */
1330 /* COMMONS USED : */
1331 /* ---------------- */
1332 /* REFERENCES CALLED : */
1333 /* ---------------------- */
1335 /* MAERMSG MCRFILL MVCVIN2 */
1338 /* DESCRIPTION/NOTES/LIMITATIONS : */
1339 /* ----------------------------------- */
1340 /* ---> Algorithm used in this general case is based on the */
1341 /* following principle : */
1342 /* Let S(t) = a0 + a1*t + a2*t**2 + ... of degree NCOEFF-1, and */
1343 /* U(t) = b0 + b1*t, then the coeff. of */
1344 /* S(U(t)) are calculated step by step with help of table TBAUX. */
1345 /* At each step number N (N=2 to NCOEFF), TBAUX(n) contains */
1346 /* the n-th coefficient of U(t)**N for n=1 to N. (RBD) */
1347 /* ---> Reference : KNUTH, 'The Art of Computer Programming', */
1348 /* Vol. 2/'Seminumerical Algorithms', */
1349 /* Ex. 11 p:451 et solution p:562. (RBD) */
1351 /* ---> Removal of the input argument CRVOLD by CRVNEW is */
1352 /* possible, which means that the call : */
1353 /* CALL MMARC41(NDIMAX,NDIMEN,NCOEFF,CURVE,UPARA0,UPARA1 */
1354 /* ,CURVE,IERCOD) */
1355 /* is absolutely LEGAL. (RBD) */
1358 /* **********************************************************************
1361 /* Name of the routine */
1363 /* Auxiliary table of coefficients of (UPARA1-UPARA0)T+UPARA0 */
1364 /* with power N=1 to NCOEFF-1. */
1367 /* Parameter adjustments */
1368 crvnew_dim1 = *ndimax;
1369 crvnew_offset = crvnew_dim1 + 1;
1370 crvnew -= crvnew_offset;
1371 crvold_dim1 = *ndimax;
1372 crvold_offset = crvold_dim1 + 1;
1373 crvold -= crvold_offset;
1377 /* **********************************************************************
1379 /* CASE WHEN PROCESSING CAN'T BE DONE */
1380 /* **********************************************************************
1382 if (*ncoeff > 61 || *ncoeff < 1) {
1386 /* **********************************************************************
1389 /* **********************************************************************
1391 if (*ndimen == *ndimax && *upara0 == 0. && *upara1 == 1.) {
1392 nboct = (*ndimax << 3) * *ncoeff;
1393 AdvApp2Var_SysBase::mcrfill_((integer *)&nboct,
1394 (char *)&crvold[crvold_offset],
1395 (char *)&crvnew[crvnew_offset]);
1398 /* **********************************************************************
1400 /* INVERSION 3D : FAST PROCESSING */
1401 /* **********************************************************************
1403 if (*upara0 == 1. && *upara1 == 0.) {
1404 if (*ndimen == 3 && *ndimax == 3 && *ncoeff <= 21) {
1405 mvcvinv_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
1409 /* ******************************************************************
1411 /* INVERSION 2D : FAST PROCESSING */
1412 /* ******************************************************************
1414 if (*ndimen == 2 && *ndimax == 2 && *ncoeff <= 21) {
1415 mvcvin2_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
1420 /* **********************************************************************
1422 /* GENERAL PROCESSING */
1423 /* **********************************************************************
1425 /* -------------------------- Initializations ---------------------------
1429 for (nd = 1; nd <= i__1; ++nd) {
1430 crvnew[nd + crvnew_dim1] = crvold[nd + crvold_dim1];
1437 tbaux[1] = *upara1 - *upara0;
1439 /* ----------------------- Calculation of coeff. of CRVNEW ------------------
1443 for (ncf = 2; ncf <= i__1; ++ncf) {
1445 /* ------------ Take into account NCF-th coeff. of CRVOLD --------
1449 for (ncj = 1; ncj <= i__2; ++ncj) {
1450 bid = tbaux[ncj - 1];
1452 for (nd = 1; nd <= i__3; ++nd) {
1453 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + ncf *
1460 bid = tbaux[ncf - 1];
1462 for (nd = 1; nd <= i__2; ++nd) {
1463 crvnew[nd + ncf * crvnew_dim1] = crvold[nd + ncf * crvold_dim1] *
1468 /* --------- Calculate (NCF+1) coeff. of ((U1-U0)*t + U0)**(NCF) ---
1471 bid = *upara1 - *upara0;
1472 tbaux[ncf] = tbaux[ncf - 1] * bid;
1473 for (ncj = ncf; ncj >= 2; --ncj) {
1474 tbaux[ncj - 1] = tbaux[ncj - 1] * *upara0 + tbaux[ncj - 2] * bid;
1477 tbaux[0] *= *upara0;
1482 /* -------------- Take into account the last coeff. of CRVOLD -----------
1486 for (ncj = 1; ncj <= i__1; ++ncj) {
1487 bid = tbaux[ncj - 1];
1489 for (nd = 1; nd <= i__2; ++nd) {
1490 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + *ncoeff *
1497 for (nd = 1; nd <= i__1; ++nd) {
1498 crvnew[nd + *ncoeff * crvnew_dim1] = crvold[nd + *ncoeff *
1499 crvold_dim1] * tbaux[*ncoeff - 1];
1503 /* ---------------------------- The end ---------------------------------
1508 AdvApp2Var_SysBase::maermsg_("MMARC41", iercod, 7L);
1514 //=======================================================================
1515 //function : AdvApp2Var_MathBase::mmarcin_
1517 //=======================================================================
1518 int AdvApp2Var_MathBase::mmarcin_(integer *ndimax,
1528 /* System generated locals */
1529 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
1533 /* Local variables */
1534 static doublereal x0, x1;
1536 static doublereal tabaux[61];
1538 static doublereal bid;
1541 static doublereal eps3;
1545 /* **********************************************************************
1548 /* Creation of curve C2(v) defined on [U0,U1] identic to */
1549 /* curve C1(u) defined on [-1,1] (change of parameter */
1550 /* of a curve) with INVERSION of indices of the resulting table. */
1554 /* GENERALIZED LIMITATION, RESTRICTION, INVERSION, CURVE */
1556 /* INPUT ARGUMENTS : */
1557 /* ------------------ */
1558 /* NDIMAX : Maximum Space Dimensioning. */
1559 /* NDIMEN : Curve Dimension. */
1560 /* NCOEFF : Nb of coefficients of the curve. */
1561 /* CRVOLD : Curve to be limited. */
1562 /* U0 : Min limit of the interval limiting the curve.
1564 /* U1 : Max limit of the interval limiting the curve.
1567 /* OUTPUT ARGUMENTS : */
1568 /* ------------------- */
1569 /* CRVNEW : Relimited curve, defined on [U0,U1] and equal to */
1570 /* CRVOLD defined on [-1,1]. */
1571 /* IERCOD : = 0, OK */
1572 /* =10, Nb of coeff. <1 or > 61. */
1573 /* =13, the requested interval of variation is null. */
1574 /* COMMONS USED : */
1575 /* ---------------- */
1576 /* REFERENCES CALLED : */
1577 /* ---------------------- */
1578 /* DESCRIPTION/NOTES/LIMITATIONS : */
1579 /* ----------------------------------- */
1581 /* **********************************************************************
1584 /* Name of the routine */
1586 /* Auxiliary table of coefficients of X1*T+X0 */
1587 /* with power N=1 to NCOEFF-1. */
1590 /* Parameter adjustments */
1591 crvnew_dim1 = *ndimax;
1592 crvnew_offset = crvnew_dim1 + 1;
1593 crvnew -= crvnew_offset;
1594 crvold_dim1 = *ncoeff;
1595 crvold_offset = crvold_dim1 + 1;
1596 crvold -= crvold_offset;
1599 ibb = AdvApp2Var_SysBase::mnfndeb_();
1601 AdvApp2Var_SysBase::mgenmsg_("MMARCIN", 7L);
1604 /* At zero machine it is tested if the output interval is not null */
1606 AdvApp2Var_MathBase::mmveps3_(&eps3);
1607 if ((d__1 = *u1 - *u0, abs(d__1)) < eps3) {
1613 /* **********************************************************************
1615 /* CASE WHEN THE PROCESSING IS IMPOSSIBLE */
1616 /* **********************************************************************
1618 if (*ncoeff > 61 || *ncoeff < 1) {
1622 /* **********************************************************************
1624 /* IF NO CHANGE OF THE INTERVAL OF DEFINITION */
1625 /* (ONLY INVERSION OF INDICES OF TABLE CRVOLD) */
1626 /* **********************************************************************
1628 if (*ndim == *ndimax && *u0 == -1. && *u1 == 1.) {
1629 AdvApp2Var_MathBase::mmcvinv_(ndim, ncoeff, ndim, &crvold[crvold_offset], &crvnew[
1633 /* **********************************************************************
1635 /* CASE WHEN THE NEW INTERVAL OF DEFINITION IS [0,1] */
1636 /* **********************************************************************
1638 if (*u0 == 0. && *u1 == 1.) {
1639 mmcvstd_(ncoeff, ndimax, ncoeff, ndim, &crvold[crvold_offset], &
1640 crvnew[crvnew_offset]);
1643 /* **********************************************************************
1645 /* GENERAL PROCESSING */
1646 /* **********************************************************************
1648 /* -------------------------- Initialization ---------------------------
1651 x0 = -(*u1 + *u0) / (*u1 - *u0);
1652 x1 = 2. / (*u1 - *u0);
1654 for (nd = 1; nd <= i__1; ++nd) {
1655 crvnew[nd + crvnew_dim1] = crvold[nd * crvold_dim1 + 1];
1664 /* ----------------------- Calculation of coeff. of CRVNEW ------------------
1668 for (ncf = 2; ncf <= i__1; ++ncf) {
1670 /* ------------ Take into account the NCF-th coeff. of CRVOLD --------
1674 for (ncj = 1; ncj <= i__2; ++ncj) {
1675 bid = tabaux[ncj - 1];
1677 for (nd = 1; nd <= i__3; ++nd) {
1678 crvnew[nd + ncj * crvnew_dim1] += crvold[ncf + nd *
1685 bid = tabaux[ncf - 1];
1687 for (nd = 1; nd <= i__2; ++nd) {
1688 crvnew[nd + ncf * crvnew_dim1] = crvold[ncf + nd * crvold_dim1] *
1693 /* --------- Calculation of (NCF+1) coeff. of [X1*t + X0]**(NCF) --------
1696 tabaux[ncf] = tabaux[ncf - 1] * x1;
1697 for (ncj = ncf; ncj >= 2; --ncj) {
1698 tabaux[ncj - 1] = tabaux[ncj - 1] * x0 + tabaux[ncj - 2] * x1;
1706 /* -------------- Take into account the last coeff. of CRVOLD -----------
1710 for (ncj = 1; ncj <= i__1; ++ncj) {
1711 bid = tabaux[ncj - 1];
1713 for (nd = 1; nd <= i__2; ++nd) {
1714 crvnew[nd + ncj * crvnew_dim1] += crvold[*ncoeff + nd *
1721 for (nd = 1; nd <= i__1; ++nd) {
1722 crvnew[nd + *ncoeff * crvnew_dim1] = crvold[*ncoeff + nd *
1723 crvold_dim1] * tabaux[*ncoeff - 1];
1727 /* ---------------------------- The end ---------------------------------
1732 AdvApp2Var_SysBase::maermsg_("MMARCIN", iercod, 7L);
1735 AdvApp2Var_SysBase::mgsomsg_("MMARCIN", 7L);
1740 //=======================================================================
1741 //function : mmatvec_
1743 //=======================================================================
1744 int mmatvec_(integer *nligne,
1755 /* System generated locals */
1758 /* Local variables */
1759 static logical ldbg;
1760 static integer jmin, jmax, i__, j, k;
1761 static doublereal somme;
1765 /* ***********************************************************************
1770 /* Produce vector matrix in form of profile */
1775 /* RESERVE, MATRIX, PRODUCT, VECTOR, PROFILE */
1777 /* INPUT ARGUMENTS : */
1778 /* -------------------- */
1779 /* NLIGNE : Line number of the matrix of constraints */
1780 /* NCOLON : Number of column of the matrix of constraints */
1781 /* GNSTOC: Number of coefficients in the profile of matrix GMATRI */
1783 /* GPOSIT: Table of positioning of terms of storage */
1784 /* GPOSIT(1,I) contains the number of terms-1 on the line I
1785 /* in the profile of the matrix. */
1786 /* GPOSIT(2,I) contains the index of storage of diagonal term*/
1788 /* GPOSIT(3,I) contains the index of column of the first term of */
1789 /* profile of line I */
1790 /* GNSTOC: Number of coefficients in the profile of matrix */
1792 /* GMATRI : Matrix of constraints in form of profile */
1793 /* VECIN : Input vector */
1794 /* DEBLIG : Line indexusing which the vector matrix is calculated */
1796 /* OUTPUT ARGUMENTS */
1797 /* --------------------- */
1798 /* VECOUT : VECTOR PRODUCT */
1800 /* IERCOD : ERROR CODE */
1803 /* COMMONS USED : */
1804 /* ------------------ */
1807 /* REFERENCES CALLED : */
1808 /* --------------------- */
1811 /* DESCRIPTION/NOTES/LIMITATIONS : */
1812 /* ----------------------------------- */
1814 /* ***********************************************************************
1817 /* ***********************************************************************
1822 /* ***********************************************************************
1824 /* INITIALISATIONS */
1825 /* ***********************************************************************
1828 /* Parameter adjustments */
1835 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1837 AdvApp2Var_SysBase::mgenmsg_("MMATVEC", 7L);
1841 /* ***********************************************************************
1844 /* ***********************************************************************
1846 AdvApp2Var_SysBase::mvriraz_((integer *)nligne,
1847 (char *)&vecout[1]);
1849 for (i__ = *deblig; i__ <= i__1; ++i__) {
1851 jmin = gposit[i__ * 3 + 3];
1852 jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
1853 aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
1855 for (j = jmin; j <= i__2; ++j) {
1857 somme += gmatri[k] * vecin[j];
1859 vecout[i__] = somme;
1868 /* ***********************************************************************
1870 /* ERROR PROCESSING */
1871 /* ***********************************************************************
1877 /* ***********************************************************************
1879 /* RETURN CALLING PROGRAM */
1880 /* ***********************************************************************
1885 /* ___ DESALLOCATION, ... */
1887 AdvApp2Var_SysBase::maermsg_("MMATVEC", iercod, 7L);
1889 AdvApp2Var_SysBase::mgsomsg_("MMATVEC", 7L);
1895 //=======================================================================
1896 //function : mmbulld_
1898 //=======================================================================
1899 int AdvApp2Var_MathBase::mmbulld_(integer *nbcoln,
1905 /* System generated locals */
1906 integer dtabtr_dim1, dtabtr_offset, i__1, i__2;
1908 /* Local variables */
1909 static logical ldbg;
1910 static doublereal daux;
1911 static integer nite1, nite2, nchan, i1, i2;
1913 /* ***********************************************************************
1918 /* Parsing of columns of a table of integers in increasing order */
1921 /* POINT-ENTRY, PARSING */
1922 /* INPUT ARGUMENTS : */
1923 /* -------------------- */
1924 /* - NBCOLN : Number of columns in the table */
1925 /* - NBLIGN : Number of lines in the table */
1926 /* - DTABTR : Table of integers to be parsed */
1927 /* - NUMCLE : Position of the key on the column */
1929 /* OUTPUT ARGUMENTS : */
1930 /* --------------------- */
1931 /* - DTABTR : Parsed table */
1933 /* COMMONS USED : */
1934 /* ------------------ */
1937 /* REFERENCES CALLED : */
1938 /* --------------------- */
1941 /* DESCRIPTION/NOTES/LIMITATIONS : */
1942 /* ----------------------------------- */
1943 /* Particularly performant if the table is almost parsed */
1944 /* In the opposite case it is better to use MVSHELD */
1945 /* ***********************************************************************
1948 /* Parameter adjustments */
1949 dtabtr_dim1 = *nblign;
1950 dtabtr_offset = dtabtr_dim1 + 1;
1951 dtabtr -= dtabtr_offset;
1954 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1956 AdvApp2Var_SysBase::mgenmsg_("MMBULLD", 7L);
1962 /* ***********************************************************************
1965 /* ***********************************************************************
1968 /* ---->ALGORITHM in N^2 / 2 additional iteration */
1972 /* ----> Parsing from left to the right */
1976 for (i1 = nite2; i1 <= i__1; ++i1) {
1977 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1 - 1)
1980 for (i2 = 1; i2 <= i__2; ++i2) {
1981 daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
1982 dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
1984 dtabtr[i2 + i1 * dtabtr_dim1] = daux;
1993 /* ----> Parsing from right to the left */
1998 for (i1 = nite1; i1 >= i__1; --i1) {
1999 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1
2000 - 1) * dtabtr_dim1]) {
2002 for (i2 = 1; i2 <= i__2; ++i2) {
2003 daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
2004 dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
2006 dtabtr[i2 + i1 * dtabtr_dim1] = daux;
2020 /* ***********************************************************************
2022 /* ERROR PROCESSING */
2023 /* ***********************************************************************
2026 /* ----> No errors at calling functions, only tests and loops. */
2028 /* ***********************************************************************
2030 /* RETURN CALLING PROGRAM */
2031 /* ***********************************************************************
2037 AdvApp2Var_SysBase::mgsomsg_("MMBULLD", 7L);
2044 //=======================================================================
2045 //function : AdvApp2Var_MathBase::mmcdriv_
2047 //=======================================================================
2048 int AdvApp2Var_MathBase::mmcdriv_(integer *ndimen,
2057 /* System generated locals */
2058 integer courbe_dim1, courbe_offset, crvdrv_dim1, crvdrv_offset, i__1,
2061 /* Local variables */
2062 static integer i__, j, k;
2063 static doublereal mfactk, bid;
2066 /* ***********************************************************************
2071 /* Calculate matrix of a derivate curve of order IDERIV. */
2072 /* with input parameters other than output parameters. */
2077 /* COEFFICIENTS,CURVE,DERIVATE I-EME. */
2079 /* INPUT ARGUMENTS : */
2080 /* ------------------ */
2081 /* NDIMEN : Space dimension (2 or 3 in general) */
2082 /* NCOEFF : Degree +1 of the curve. */
2083 /* COURBE : Table of coefficients of the curve. */
2084 /* IDERIV : Required order of derivation : 1=1st derivate, etc... */
2086 /* OUTPUT ARGUMENTS : */
2087 /* ------------------- */
2088 /* NCOFDV : Degree +1 of the derivative of order IDERIV of the curve. */
2089 /* CRVDRV : Table of coefficients of the derivative of order IDERIV */
2092 /* COMMONS USED : */
2093 /* ---------------- */
2095 /* REFERENCES CALLED : */
2096 /* ----------------------- */
2098 /* DESCRIPTION/NOTES/LIMITATIONS : */
2099 /* ----------------------------------- */
2101 /* ---> It is possible to take as output argument the curve */
2102 /* and the number of coeff passed at input by making : */
2103 /* CALL MMCDRIV(NDIMEN,NCOEFF,COURBE,IDERIV,NCOEFF,COURBE). */
2104 /* After this call, NCOEFF does the number of coeff of the derived */
2105 /* curve the coefficients which of are stored in CURVE. */
2106 /* Attention to the coefficients of CURVE of rank superior to */
2107 /* NCOEFF : they are not set to zero. */
2109 /* ---> Algorithm : */
2110 /* The code below was written basing on the following algorithm:
2113 /* Let P(t) = a1 + a2*t + ... an*t**n. Derivate of order k of P */
2114 /* (containing n-k coefficients) is calculated as follows : */
2116 /* Pk(t) = a(k+1)*CNP(k,k)*k! */
2117 /* + a(k+2)*CNP(k+1,k)*k! * t */
2121 /* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
2122 /* ***********************************************************************
2126 /* -------------- Case when the order of derivative is -------------------
2128 /* ---------------- greater than the degree of the curve ---------------------
2131 /* **********************************************************************
2136 /* Serves to provide the coefficients of binome (Pascal's triangle). */
2140 /* Binomial coeff from 0 to 60. read only . init par block data */
2142 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
2143 /* ----------------------------------- */
2144 /* Binomial coefficients form a triangular matrix. */
2145 /* This matrix is completed in table CNP by its transposition. */
2146 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
2148 /* Initialization is done by block-data MMLLL09.RES, */
2149 /* created by program MQINICNP.FOR). */
2150 /* **********************************************************************
2155 /* ***********************************************************************
2158 /* Parameter adjustments */
2159 crvdrv_dim1 = *ndimen;
2160 crvdrv_offset = crvdrv_dim1 + 1;
2161 crvdrv -= crvdrv_offset;
2162 courbe_dim1 = *ndimen;
2163 courbe_offset = courbe_dim1 + 1;
2164 courbe -= courbe_offset;
2167 if (*ideriv >= *ncoeff) {
2169 for (i__ = 1; i__ <= i__1; ++i__) {
2170 crvdrv[i__ + crvdrv_dim1] = 0.;
2176 /* **********************************************************************
2178 /* General processing */
2179 /* **********************************************************************
2181 /* --------------------- Calculation of Factorial(IDERIV) ------------------
2187 for (i__ = 2; i__ <= i__1; ++i__) {
2192 /* ------------ Calculation of coeff of the derived of order IDERIV ----------
2194 /* ---> Attention : coefficient binomial C(n,m) is represented in */
2195 /* MCCNP by CNP(N+1,M+1). */
2198 for (j = k + 1; j <= i__1; ++j) {
2199 bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
2201 for (i__ = 1; i__ <= i__2; ++i__) {
2202 crvdrv[i__ + (j - k) * crvdrv_dim1] = bid * courbe[i__ + j *
2209 *ncofdv = *ncoeff - *ideriv;
2211 /* -------------------------------- The end -----------------------------
2218 //=======================================================================
2219 //function : AdvApp2Var_MathBase::mmcglc1_
2221 //=======================================================================
2222 int AdvApp2Var_MathBase::mmcglc1_(integer *ndimax,
2235 /* System generated locals */
2236 integer courbe_dim1, courbe_offset, i__1;
2239 /* Local variables */
2240 static integer ndec;
2241 static doublereal tdeb, tfin;
2242 static integer iter;
2243 static doublereal oldso;
2244 static integer itmax;
2245 static doublereal sottc;
2246 static integer kk, ibb;
2247 static doublereal dif, pas;
2248 static doublereal som;
2251 /* ***********************************************************************
2256 /* Allows calculating the length of an arc of curve POLYNOMIAL */
2257 /* on an interval [A,B]. */
2261 /* LENGTH,CURVE,GAUSS,PRIVATE. */
2263 /* INPUT ARGUMENTS : */
2264 /* ------------------ */
2265 /* NDIMAX : Max. number of lines of tables */
2266 /* (i.e. max. nb of polynoms). */
2267 /* NDIMEN : Dimension of the space (nb of polynoms). */
2268 /* NCOEFF : Nb of coefficients of the polynom. This is degree + 1.
2270 /* COURBE(NDIMAX,NCOEFF) : Coefficients of the curve. */
2271 /* TDEBUT : Lower limit of the interval of integration for */
2272 /* length calculation. */
2273 /* TFINAL : Upper limit of the interval of integration for */
2274 /* length calculation. */
2275 /* EPSILN : REQIRED precision for length calculation. */
2277 /* OUTPUT ARGUMENTS : */
2278 /* ------------------- */
2279 /* XLONGC : Length of the arc of curve */
2280 /* ERREUR : Precision OBTAINED for the length calculation. */
2281 /* IERCOD : Error code, 0 OK, >0 Serious error. */
2282 /* = 1 Too much iterations, the best calculated resultat */
2283 /* (is almost ERROR) */
2284 /* = 2 Pb MMLONCV (no result) */
2285 /* = 3 NDIM or NCOEFF invalid (no result) */
2287 /* COMMONS USED : */
2288 /* ---------------- */
2290 /* REFERENCES CALLED : */
2291 /* ----------------------- */
2293 /* DESCRIPTION/NOTES/LIMITATIONS : */
2294 /* ----------------------------------- */
2295 /* The polynom is actually a set of polynoms with */
2296 /* coefficients arranged in a table of 2 indices, */
2297 /* each line relative to the polynom. */
2298 /* The polynom is defined by these coefficients ordered */
2299 /* by increasing power of the variable. */
2300 /* All polynoms have the same number of coefficients (the */
2303 /* This program cancels and replaces LENGCV, MLONGC and MLENCV. */
2305 /* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
2308 /* ***********************************************************************
2311 /* Name of the routine */
2314 /* ------------------------ General Initialization ---------------------
2317 /* Parameter adjustments */
2318 courbe_dim1 = *ndimax;
2319 courbe_offset = courbe_dim1 + 1;
2320 courbe -= courbe_offset;
2323 ibb = AdvApp2Var_SysBase::mnfndeb_();
2325 AdvApp2Var_SysBase::mgenmsg_("MMCGLC1", 7L);
2332 /* ------ Test of equity of limits */
2334 if (*tdebut == *tfinal) {
2339 /* ------ Test of the dimension and the number of coefficients */
2341 if (*ndimen <= 0 || *ncoeff <= 0) {
2345 /* ----- Nb of current cutting, nb of iteration, */
2346 /* max nb of iterations */
2353 /* ------ Variation of the nb of intervals */
2354 /* Multiplied by 2 at each iteration */
2357 pas = (*tfinal - *tdebut) / ndec;
2360 /* ------ Loop on all current NDEC intervals */
2363 for (kk = 1; kk <= i__1; ++kk) {
2365 /* ------ Limits of the current integration interval */
2367 tdeb = *tdebut + (kk - 1) * pas;
2369 mmloncv_(ndimax, ndimen, ncoeff, &courbe[courbe_offset], &tdeb, &tfin,
2381 /* ----------------- Test of the maximum number of iterations ------------
2384 /* Test if passes at least once ** */
2393 /* ------ Take into account DIF - Test of convergence */
2396 dif = (d__1 = sottc - oldso, abs(d__1));
2398 /* ------ If DIF is OK, leave..., otherwise: */
2400 if (dif > *epsiln) {
2402 /* ------ If nb iteration exceeded, leave */
2409 /* ------ Otherwise continue by cutting the initial interval.
2419 /* ------------------------------ THE END -------------------------------
2427 /* ---> PB in MMLONCV */
2433 /* ---> NCOEFF or NDIM invalid. */
2441 AdvApp2Var_SysBase::maermsg_("MMCGLC1", iercod, 7L);
2444 AdvApp2Var_SysBase::mgsomsg_("MMCGLC1", 7L);
2449 //=======================================================================
2450 //function : mmchole_
2452 //=======================================================================
2453 int mmchole_(integer *,//mxcoef,
2462 /* System generated locals */
2463 integer i__1, i__2, i__3;
2466 /* Builtin functions */
2469 /* Local variables */
2470 static logical ldbg;
2471 static integer kmin, i__, j, k;
2472 static doublereal somme;
2473 static integer ptini, ptcou;
2476 /* ***********************************************************************
2481 /* Produce decomposition of choleski of matrix A in S.S */
2482 /* Calculate inferior triangular matrix S. */
2486 /* RESOLUTION, MFACTORISATION, MATRIX_PROFILE, CHOLESKI */
2488 /* INPUT ARGUMENTS : */
2489 /* -------------------- */
2490 /* MXCOEF : Max number of terms in the hessian profile */
2491 /* DIMENS : Dimension of the problem */
2492 /* AMATRI(MXCOEF) : Coefficients of the matrix profile */
2493 /* APOSIT(1,*) : Distance diagonal-left extremity of the line
2495 /* APOSIT(2,*) : Position of diagonal terms in HESSIE */
2496 /* POSUIV(MXCOEF) : first line inferior not out of profile */
2498 /* OUTPUT ARGUMENTS : */
2499 /* --------------------- */
2500 /* CHOMAT(MXCOEF) : Inferior triangular matrix preserving the */
2501 /* profile of AMATRI. */
2502 /* IERCOD : error code */
2504 /* = 1 : non-defined positive matrix */
2506 /* COMMONS USED : */
2507 /* ------------------ */
2511 /* REFERENCES CALLED : */
2512 /* ---------------------- */
2514 /* DESCRIPTION/NOTES/LIMITATIONS : */
2515 /* ----------------------------------- */
2516 /* DEBUG LEVEL = 4 */
2517 /* ***********************************************************************
2520 /* ***********************************************************************
2525 /* ***********************************************************************
2527 /* INITIALISATIONS */
2528 /* ***********************************************************************
2531 /* Parameter adjustments */
2538 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
2540 AdvApp2Var_SysBase::mgenmsg_("MMCHOLE", 7L);
2544 /* ***********************************************************************
2547 /* ***********************************************************************
2551 for (j = 1; j <= i__1; ++j) {
2553 ptini = aposit[(j << 1) + 2];
2557 for (k = ptini - aposit[(j << 1) + 1]; k <= i__2; ++k) {
2558 /* Computing 2nd power */
2560 somme += d__1 * d__1;
2563 if (amatri[ptini] - somme < 1e-32) {
2566 chomat[ptini] = sqrt(amatri[ptini] - somme);
2570 while(posuiv[ptcou] > 0) {
2572 i__ = posuiv[ptcou];
2573 ptcou = aposit[(i__ << 1) + 2] - (i__ - j);
2575 /* Calculate the sum of S .S for k =1 a j-1 */
2579 i__2 = i__ - aposit[(i__ << 1) + 1], i__3 = j - aposit[(j << 1) +
2581 kmin = max(i__2,i__3);
2583 for (k = kmin; k <= i__2; ++k) {
2584 somme += chomat[aposit[(i__ << 1) + 2] - (i__ - k)] * chomat[
2585 aposit[(j << 1) + 2] - (j - k)];
2588 chomat[ptcou] = (amatri[ptcou] - somme) / chomat[ptini];
2594 /* ***********************************************************************
2596 /* ERROR PROCESSING */
2597 /* ***********************************************************************
2604 /* ***********************************************************************
2606 /* RETURN CALLING PROGRAM */
2607 /* ***********************************************************************
2612 AdvApp2Var_SysBase::maermsg_("MMCHOLE", iercod, 7L);
2614 AdvApp2Var_SysBase::mgsomsg_("MMCHOLE", 7L);
2620 //=======================================================================
2621 //function : AdvApp2Var_MathBase::mmcvctx_
2623 //=======================================================================
2624 int AdvApp2Var_MathBase::mmcvctx_(integer *ndimen,
2634 /* System generated locals */
2635 integer ctrtes_dim1, ctrtes_offset, crvres_dim1, crvres_offset,
2636 xmatri_dim1, xmatri_offset, tabaux_dim1, tabaux_offset, i__1,
2639 /* Local variables */
2640 static integer moup1, nordr;
2642 static integer ibb, ncf, ndv;
2643 static doublereal eps1;
2646 /* ***********************************************************************
2651 /* Calculate a polynomial curve checking the */
2652 /* passage constraints (interpolation) */
2653 /* from first derivatives, etc... to extremities. */
2654 /* Parameters at the extremities are supposed to be -1 and 1. */
2658 /* ALL, AB_SPECIFI::CONSTRAINTS&,INTERPOLATION,&CURVE */
2660 /* INPUT ARGUMENTS : */
2661 /* ------------------ */
2662 /* NDIMEN : Space Dimension. */
2663 /* NCOFMX : Nb of coeff. of curve CRVRES on each */
2665 /* NDERIV : Order of constraint with derivatives : */
2666 /* 0 --> interpolation simple. */
2667 /* 1 --> interpolation+constraints with 1st. */
2668 /* 2 --> cas (0)+ (1) + " " 2nd derivatives. */
2670 /* CTRTES : Table of constraints. */
2671 /* CTRTES(*,1,*) = contraints at -1. */
2672 /* CTRTES(*,2,*) = contraints at 1. */
2674 /* OUTPUT ARGUMENTS : */
2675 /* ------------------- */
2676 /* CRVRES : Resulting curve defined on (-1,1). */
2677 /* TABAUX : Auxilliary matrix. */
2678 /* XMATRI : Auxilliary matrix. */
2680 /* COMMONS UTILISES : */
2681 /* ---------------- */
2685 /* REFERENCES CALLED : */
2686 /* ---------------------- */
2688 /* MAERMSG R*8 DFLOAT MGENMSG */
2689 /* MGSOMSG MMEPS1 MMRSLW */
2692 /* DESCRIPTION/NOTES/LIMITATIONS : */
2693 /* ----------------------------------- */
2694 /* The polynom (or the curve) is calculated by solving a */
2695 /* system of linear equations. If the imposed degree is great */
2696 /* it is preferable to call a routine based on */
2697 /* Lagrange or Hermite interpolation depending on the case. */
2698 /* (for a high degree the matrix of the system can be badly */
2699 /* conditionned). */
2700 /* This routine returns a curve defined in (-1,1). */
2701 /* In general case, it is necessary to use MCVCTG. */
2703 /* ***********************************************************************
2706 /* Name of the routine */
2709 /* Parameter adjustments */
2710 crvres_dim1 = *ncofmx;
2711 crvres_offset = crvres_dim1 + 1;
2712 crvres -= crvres_offset;
2713 xmatri_dim1 = *nderiv + 1;
2714 xmatri_offset = xmatri_dim1 + 1;
2715 xmatri -= xmatri_offset;
2716 tabaux_dim1 = *nderiv + 1 + *ndimen;
2717 tabaux_offset = tabaux_dim1 + 1;
2718 tabaux -= tabaux_offset;
2719 ctrtes_dim1 = *ndimen;
2720 ctrtes_offset = ctrtes_dim1 * 3 + 1;
2721 ctrtes -= ctrtes_offset;
2724 ibb = AdvApp2Var_SysBase::mnfndeb_();
2726 AdvApp2Var_SysBase::mgenmsg_("MMCVCTX", 7L);
2729 AdvApp2Var_MathBase::mmeps1_(&eps1);
2731 /* ****************** CALCULATION OF EVEN COEFFICIENTS *********************
2733 /* ------------------------- Initialization -----------------------------
2736 nordr = *nderiv + 1;
2738 for (ncf = 1; ncf <= i__1; ++ncf) {
2739 tabaux[ncf + tabaux_dim1] = 1.;
2743 /* ---------------- Calculation of terms corresponding to derivatives -------
2747 for (ndv = 2; ndv <= i__1; ++ndv) {
2749 for (ncf = 1; ncf <= i__2; ++ncf) {
2750 tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
2751 tabaux_dim1] * (doublereal) ((ncf << 1) - ndv);
2757 /* ------------------ Writing the second member -----------------------
2762 for (ndv = 1; ndv <= i__1; ++ndv) {
2764 for (nd = 1; nd <= i__2; ++nd) {
2765 tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
2766 + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2767 * ctrtes_dim1]) / 2.;
2774 /* -------------------- Resolution of the system ---------------------------
2777 mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2778 xmatri_offset], iercod);
2783 for (nd = 1; nd <= i__1; ++nd) {
2785 for (ncf = 1; ncf <= i__2; ++ncf) {
2786 crvres[(ncf << 1) - 1 + nd * crvres_dim1] = xmatri[ncf + nd *
2793 /* ***************** CALCULATION OF UNEVEN COEFFICIENTS ********************
2795 /* ------------------------- Initialization -----------------------------
2800 for (ncf = 1; ncf <= i__1; ++ncf) {
2801 tabaux[ncf + tabaux_dim1] = 1.;
2805 /* ---------------- Calculation of terms corresponding to derivatives -------
2809 for (ndv = 2; ndv <= i__1; ++ndv) {
2811 for (ncf = 1; ncf <= i__2; ++ncf) {
2812 tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
2813 tabaux_dim1] * (doublereal) ((ncf << 1) - ndv + 1);
2819 /* ------------------ Writing of the second member -----------------------
2824 for (ndv = 1; ndv <= i__1; ++ndv) {
2826 for (nd = 1; nd <= i__2; ++nd) {
2827 tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
2828 + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2829 * ctrtes_dim1]) / 2.;
2836 /* -------------------- Solution of the system ---------------------------
2839 mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2840 xmatri_offset], iercod);
2845 for (nd = 1; nd <= i__1; ++nd) {
2847 for (ncf = 1; ncf <= i__2; ++ncf) {
2848 crvres[(ncf << 1) + nd * crvres_dim1] = xmatri[ncf + nd *
2855 /* --------------------------- The end ----------------------------------
2860 AdvApp2Var_SysBase::maermsg_("MMCVCTX", iercod, 7L);
2863 AdvApp2Var_SysBase::mgsomsg_("MMCVCTX", 7L);
2869 //=======================================================================
2870 //function : AdvApp2Var_MathBase::mmcvinv_
2872 //=======================================================================
2873 int AdvApp2Var_MathBase::mmcvinv_(integer *ndimax,
2880 /* Initialized data */
2882 static char nomprg[8+1] = "MMCVINV ";
2884 /* System generated locals */
2885 integer curve_dim1, curve_offset, curveo_dim1, curveo_offset, i__1, i__2;
2887 /* Local variables */
2888 static integer i__, nd, ibb;
2891 /* ***********************************************************************
2896 /* Inversion of arguments of the final curve. */
2900 /* SMOOTHING,CURVE */
2903 /* INPUT ARGUMENTS : */
2904 /* ------------------ */
2906 /* NDIM: Space Dimension. */
2907 /* NCOEF: Degree of the polynom. */
2908 /* CURVEO: The curve before inversion. */
2910 /* OUTPUT ARGUMENTS : */
2911 /* ------------------- */
2912 /* CURVE: The curve after inversion. */
2914 /* COMMONS USED : */
2915 /* ---------------- */
2916 /* REFERENCES APPELEES : */
2917 /* ----------------------- */
2918 /* DESCRIPTION/NOTES/LIMITATIONS : */
2919 /* ----------------------------------- */
2920 /* ***********************************************************************
2923 /* The name of the routine */
2924 /* Parameter adjustments */
2925 curve_dim1 = *ndimax;
2926 curve_offset = curve_dim1 + 1;
2927 curve -= curve_offset;
2928 curveo_dim1 = *ncoef;
2929 curveo_offset = curveo_dim1 + 1;
2930 curveo -= curveo_offset;
2934 ibb = AdvApp2Var_SysBase::mnfndeb_();
2936 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
2940 for (i__ = 1; i__ <= i__1; ++i__) {
2942 for (nd = 1; nd <= i__2; ++nd) {
2943 curve[nd + i__ * curve_dim1] = curveo[i__ + nd * curveo_dim1];
2952 //=======================================================================
2953 //function : mmcvstd_
2955 //=======================================================================
2956 int mmcvstd_(integer *ncofmx,
2964 /* System generated locals */
2965 integer courbe_dim1, crvcan_dim1, crvcan_offset, i__1, i__2, i__3;
2967 /* Local variables */
2968 static integer ndeg, i__, j, j1, nd, ibb;
2969 static doublereal bid;
2972 /* ***********************************************************************
2977 /* Transform curve defined between [-1,1] into [0,1]. */
2981 /* LIMITATION,RESTRICTION,CURVE */
2983 /* INPUT ARGUMENTS : */
2984 /* ------------------ */
2985 /* NDIMAX : Dimension of the space. */
2986 /* NDIMEN : Dimension of the curve. */
2987 /* NCOEFF : Degree of the curve. */
2988 /* CRVCAN(NCOFMX,NDIMEN): The curve is defined at the interval [-1,1]. */
2990 /* OUTPUT ARGUMENTS : */
2991 /* ------------------- */
2992 /* CURVE(NDIMAX,NCOEFF): Curve defined at the interval [0,1]. */
2994 /* COMMONS USED : */
2995 /* ---------------- */
2997 /* REFERENCES CALLED : */
2998 /* ----------------------- */
3000 /* DESCRIPTION/NOTES/LIMITATIONS : */
3001 /* ----------------------------------- */
3003 /* ***********************************************************************
3006 /* Name of the program. */
3009 /* **********************************************************************
3014 /* Provides binomial coefficients (Pascal triangle). */
3018 /* Binomial coefficient from 0 to 60. read only . init by block data */
3020 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3021 /* ----------------------------------- */
3022 /* Binomial coefficients form a triangular matrix. */
3023 /* This matrix is completed in table CNP by its transposition. */
3024 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
3026 /* Initialization is done with block-data MMLLL09.RES, */
3027 /* created by the program MQINICNP.FOR. */
3029 /* **********************************************************************
3034 /* ***********************************************************************
3037 /* Parameter adjustments */
3038 courbe_dim1 = *ndimax;
3040 crvcan_dim1 = *ncofmx;
3041 crvcan_offset = crvcan_dim1;
3042 crvcan -= crvcan_offset;
3045 ibb = AdvApp2Var_SysBase::mnfndeb_();
3047 AdvApp2Var_SysBase::mgenmsg_("MMCVSTD", 7L);
3051 /* ------------------ Construction of the resulting curve ----------------
3055 for (nd = 1; nd <= i__1; ++nd) {
3057 for (j = 0; j <= i__2; ++j) {
3060 for (i__ = j; i__ <= i__3; i__ += 2) {
3061 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
3065 courbe[nd + j * courbe_dim1] = bid;
3070 for (i__ = j1; i__ <= i__3; i__ += 2) {
3071 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
3075 courbe[nd + j * courbe_dim1] -= bid;
3081 /* ------------------- Renormalization of the CURVE -------------------------
3086 for (i__ = 0; i__ <= i__1; ++i__) {
3088 for (nd = 1; nd <= i__2; ++nd) {
3089 courbe[nd + i__ * courbe_dim1] *= bid;
3096 /* ----------------------------- The end --------------------------------
3100 AdvApp2Var_SysBase::mgsomsg_("MMCVSTD", 7L);
3105 //=======================================================================
3106 //function : AdvApp2Var_MathBase::mmdrc11_
3108 //=======================================================================
3109 int AdvApp2Var_MathBase::mmdrc11_(integer *iordre,
3114 doublereal *mfactab)
3117 /* System generated locals */
3118 integer courbe_dim1, courbe_offset, points_dim2, points_offset, i__1,
3121 /* Local variables */
3123 static integer ndeg, i__, j, ndgcb, nd, ibb;
3126 /* **********************************************************************
3131 /* Calculation of successive derivatives of equation CURVE with */
3132 /* parameters -1, 1 from order 0 to order IORDRE */
3133 /* included. The calculation is produced without knowing the coefficients of */
3134 /* derivatives of the curve. */
3138 /* POSITIONING,EXTREMITIES,CURVE,DERIVATIVE. */
3140 /* INPUT ARGUMENTS : */
3141 /* ------------------ */
3142 /* IORDRE : Maximum order of calculation of derivatives. */
3143 /* NDIMEN : Dimension of the space. */
3144 /* NCOEFF : Number of coefficients of the curve (degree+1). */
3145 /* COURBE : Table of coefficients of the curve. */
3147 /* OUTPUT ARGUMENTS : */
3148 /* ------------------- */
3149 /* POINTS : Table of values of consecutive derivatives */
3150 /* of parameters -1.D0 and 1.D0. */
3151 /* MFACTAB : Auxiliary table for calculation of factorial(I).
3154 /* COMMONS USED : */
3155 /* ---------------- */
3158 /* REFERENCES CALLED : */
3159 /* ----------------------- */
3161 /* DESCRIPTION/NOTES/LIMITATIONS : */
3162 /* ----------------------------------- */
3164 /* ---> ATTENTION, the coefficients of the curve are */
3165 /* in a reverse order. */
3167 /* ---> The algorithm of calculation of derivatives is based on */
3168 /* generalization of Horner scheme : */
3170 /* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
3173 /* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
3175 /* aj = a(j-1).x + u(k-j) */
3176 /* bj = b(j-1).x + a(j-1) */
3177 /* cj = c(j-1).x + b(j-1) */
3179 /* So : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
3181 /* The algorithm is generalized easily for calculation of */
3188 /* Reference : D. KNUTH, "The Art of Computer Programming" */
3189 /* --------- Vol. 2/Seminumerical Algorithms */
3190 /* Addison-Wesley Pub. Co. (1969) */
3191 /* pages 423-425. */
3193 /* **********************************************************************
3196 /* Name of the routine */
3198 /* Parameter adjustments */
3199 points_dim2 = *iordre + 1;
3200 points_offset = (points_dim2 << 1) + 1;
3201 points -= points_offset;
3202 courbe_dim1 = *ncoeff;
3203 courbe_offset = courbe_dim1;
3204 courbe -= courbe_offset;
3207 ibb = AdvApp2Var_SysBase::mnfndeb_();
3209 AdvApp2Var_SysBase::mgenmsg_("MMDRC11", 7L);
3212 if (*iordre < 0 || *ncoeff < 1) {
3216 /* ------------------- Initialization of table POINTS -----------------
3219 ndgcb = *ncoeff - 1;
3221 for (nd = 1; nd <= i__1; ++nd) {
3222 points[(nd * points_dim2 << 1) + 1] = courbe[ndgcb + nd * courbe_dim1]
3224 points[(nd * points_dim2 << 1) + 2] = courbe[ndgcb + nd * courbe_dim1]
3230 for (nd = 1; nd <= i__1; ++nd) {
3232 for (j = 1; j <= i__2; ++j) {
3233 points[((j + nd * points_dim2) << 1) + 1] = 0.;
3234 points[((j + nd * points_dim2) << 1) + 2] = 0.;
3240 /* Calculation with parameter -1 and 1 */
3243 for (nd = 1; nd <= i__1; ++nd) {
3245 for (ndeg = 1; ndeg <= i__2; ++ndeg) {
3246 for (i__ = *iordre; i__ >= 1; --i__) {
3247 points[((i__ + nd * points_dim2) << 1) + 1] = -points[((i__ + nd
3248 * points_dim2) << 1) + 1] + points[((i__ - 1 + nd *
3249 points_dim2) << 1) + 1];
3250 points[((i__ + nd * points_dim2) << 1) + 2] += points[((i__ - 1
3251 + nd * points_dim2) << 1) + 2];
3254 points[(nd * points_dim2 << 1) + 1] = -points[(nd * points_dim2 <<
3255 1) + 1] + courbe[ndgcb - ndeg + nd * courbe_dim1];
3256 points[(nd * points_dim2 << 1) + 2] += courbe[ndgcb - ndeg + nd *
3263 /* --------------------- Multiplication by factorial(I) --------------
3267 mfac_(&mfactab[1], iordre);
3270 for (nd = 1; nd <= i__1; ++nd) {
3272 for (i__ = 2; i__ <= i__2; ++i__) {
3273 points[((i__ + nd * points_dim2) << 1) + 1] = mfactab[i__] *
3274 points[((i__ + nd * points_dim2) << 1) + 1];
3275 points[((i__ + nd * points_dim2) << 1) + 2] = mfactab[i__] *
3276 points[((i__ + nd * points_dim2) << 1) + 2];
3283 /* ---------------------------- End -------------------------------------
3288 AdvApp2Var_SysBase::mgsomsg_("MMDRC11", 7L);
3293 //=======================================================================
3294 //function : mmdrvcb_
3296 //=======================================================================
3297 int mmdrvcb_(integer *ideriv,
3306 /* System generated locals */
3307 integer courbe_dim1, tabpnt_dim1, i__1, i__2, i__3;
3309 /* Local variables */
3310 static integer ndeg, i__, j, nd, ndgcrb, iptpnt, ibb;
3313 /* ***********************************************************************
3317 /* Calculation of successive derivatives of equation CURVE with */
3318 /* parameter TPARAM from order 0 to order IDERIV included. */
3319 /* The calculation is produced without knowing the coefficients of */
3320 /* derivatives of the CURVE. */
3324 /* POSITIONING,PARAMETER,CURVE,DERIVATIVE. */
3326 /* INPUT ARGUMENTS : */
3327 /* ------------------ */
3328 /* IORDRE : Maximum order of calculation of derivatives. */
3329 /* NDIMEN : Dimension of the space. */
3330 /* NCOEFF : Number of coefficients of the curve (degree+1). */
3331 /* COURBE : Table of coefficients of the curve. */
3332 /* TPARAM : Value of the parameter where the curve should be evaluated. */
3334 /* OUTPUT ARGUMENTS : */
3335 /* ------------------- */
3336 /* TABPNT : Table of values of consecutive derivatives */
3337 /* of parameter TPARAM. */
3338 /* IERCOD : 0 = OK, */
3339 /* 1 = incoherent input. */
3341 /* COMMONS USED : */
3342 /* ---------------- */
3345 /* REFERENCES CALLED : */
3346 /* ----------------------- */
3348 /* DESCRIPTION/NOTES/LIMITATIONS : */
3349 /* ----------------------------------- */
3351 /* The algorithm of calculation of derivatives is based on */
3352 /* generalization of the Horner scheme : */
3354 /* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
3357 /* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
3359 /* aj = a(j-1).x + u(k-j) */
3360 /* bj = b(j-1).x + a(j-1) */
3361 /* cj = c(j-1).x + b(j-1) */
3363 /* So, it is obtained : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
3365 /* The algorithm can be easily generalized for the calculation of */
3372 /* Reference : D. KNUTH, "The Art of Computer Programming" */
3373 /* --------- Vol. 2/Seminumerical Algorithms */
3374 /* Addison-Wesley Pub. Co. (1969) */
3375 /* pages 423-425. */
3377 /* ---> To evaluare derivatives at 0 and 1, it is preferable */
3378 /* to use routine MDRV01.FOR . */
3380 /* **********************************************************************
3383 /* Name of the routine */
3385 /* Parameter adjustments */
3386 tabpnt_dim1 = *ndim;
3388 courbe_dim1 = *ndim;
3392 ibb = AdvApp2Var_SysBase::mnfndeb_();
3394 AdvApp2Var_SysBase::mgenmsg_("MMDRVCB", 7L);
3397 if (*ideriv < 0 || *ncoeff < 1) {
3403 /* ------------------- Initialization of table TABPNT -----------------
3406 ndgcrb = *ncoeff - 1;
3408 for (nd = 1; nd <= i__1; ++nd) {
3409 tabpnt[nd] = courbe[nd + ndgcrb * courbe_dim1];
3416 iptpnt = *ndim * *ideriv;
3417 AdvApp2Var_SysBase::mvriraz_((integer *)&iptpnt,
3418 (char *)&tabpnt[tabpnt_dim1 + 1]);
3421 /* ------------------------ Calculation of parameter TPARAM ------------------
3425 for (ndeg = 1; ndeg <= i__1; ++ndeg) {
3427 for (nd = 1; nd <= i__2; ++nd) {
3428 for (i__ = *ideriv; i__ >= 1; --i__) {
3429 tabpnt[nd + i__ * tabpnt_dim1] = tabpnt[nd + i__ *
3430 tabpnt_dim1] * *tparam + tabpnt[nd + (i__ - 1) *
3434 tabpnt[nd] = tabpnt[nd] * *tparam + courbe[nd + (ndgcrb - ndeg) *
3441 /* --------------------- Multiplication by factorial(I) -------------
3445 for (i__ = 2; i__ <= i__1; ++i__) {
3447 for (j = 2; j <= i__2; ++j) {
3449 for (nd = 1; nd <= i__3; ++nd) {
3450 tabpnt[nd + i__ * tabpnt_dim1] = (doublereal) j * tabpnt[nd +
3459 /* --------------------------- The end ---------------------------------
3464 AdvApp2Var_SysBase::maermsg_("MMDRVCB", iercod, 7L);
3469 //=======================================================================
3470 //function : AdvApp2Var_MathBase::mmdrvck_
3472 //=======================================================================
3473 int AdvApp2Var_MathBase::mmdrvck_(integer *ncoeff,
3481 /* Initialized data */
3483 static doublereal mmfack[21] = { 1.,2.,6.,24.,120.,720.,5040.,40320.,
3484 362880.,3628800.,39916800.,479001600.,6227020800.,87178291200.,
3485 1.307674368e12,2.0922789888e13,3.55687428096e14,6.402373705728e15,
3486 1.21645100408832e17,2.43290200817664e18,5.109094217170944e19 };
3488 /* System generated locals */
3489 integer courbe_dim1, courbe_offset, i__1, i__2;
3491 /* Local variables */
3492 static integer i__, j, k, nd;
3493 static doublereal mfactk, bid;
3496 /* IMPLICIT INTEGER (I-N) */
3497 /* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
3500 /* ***********************************************************************
3505 /* Calculate the value of a derived curve of order IDERIV in */
3506 /* a point of parameter TPARAM. */
3510 /* POSITIONING,CURVE,DERIVATIVE of ORDER K. */
3512 /* INPUT ARGUMENTS : */
3513 /* ------------------ */
3514 /* NCOEFF : Degree +1 of the curve. */
3515 /* NDIMEN : Dimension of the space (2 or 3 in general) */
3516 /* COURBE : Table of coefficients of the curve. */
3517 /* IDERIV : Required order of derivation : 1=1st derivative, etc... */
3518 /* TPARAM : Value of parameter of the curve. */
3520 /* OUTPUT ARGUMENTS : */
3521 /* ------------------- */
3522 /* PNTCRB : Point of parameter TPARAM on the derivative of order */
3523 /* IDERIV of CURVE. */
3525 /* COMMONS USED : */
3526 /* ---------------- */
3529 /* REFERENCES CALLED : */
3530 /* ---------------------- */
3532 /* DESCRIPTION/NOTES/LIMITATIONS : */
3533 /* ----------------------------------- */
3535 /* The code below was written basing on the following algorithm :
3538 /* Let P(t) = a1 + a2*t + ... an*t**n. The derivative of order k of P */
3539 /* (containing n-k coefficients) is calculated as follows : */
3541 /* Pk(t) = a(k+1)*CNP(k,k)*k! */
3542 /* + a(k+2)*CNP(k+1,k)*k! * t */
3546 /* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
3548 /* Evaluation is produced following the classic Horner scheme. */
3550 /* ***********************************************************************
3554 /* Factorials (1 to 21) caculated on VAX in R*16 */
3557 /* **********************************************************************
3562 /* Serves to provide binomial coefficients (Pascal triangle). */
3566 /* Binomial Coeff from 0 to 60. read only . init by block data */
3568 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3569 /* ----------------------------------- */
3570 /* Binomial coefficients form a triangular matrix. */
3571 /* This matrix is completed in table CNP by its transposition. */
3572 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
3574 /* Initialization is done by block-data MMLLL09.RES, */
3575 /* created by program MQINICNP.FOR. */
3577 /* **********************************************************************
3582 /* ***********************************************************************
3585 /* Parameter adjustments */
3587 courbe_dim1 = *ndimen;
3588 courbe_offset = courbe_dim1 + 1;
3589 courbe -= courbe_offset;
3593 /* -------------- Case when the order of derivative is greater than -------------------
3595 /* ---------------- the degree of the curve ---------------------
3598 if (*ideriv >= *ncoeff) {
3600 for (nd = 1; nd <= i__1; ++nd) {
3606 /* **********************************************************************
3608 /* General processing*/
3609 /* **********************************************************************
3611 /* --------------------- Calculation of Factorial(IDERIV) ------------------
3615 if (*ideriv <= 21 && *ideriv > 0) {
3616 mfactk = mmfack[k - 1];
3620 for (i__ = 2; i__ <= i__1; ++i__) {
3626 /* ------- Calculation of derivative of order IDERIV of CURVE in TPARAM -----
3628 /* ---> Attention : binomial coefficient C(n,m) is represented in */
3629 /* MCCNP by CNP(N,M). */
3632 for (nd = 1; nd <= i__1; ++nd) {
3633 pntcrb[nd] = courbe[nd + *ncoeff * courbe_dim1] * mmcmcnp_.cnp[*
3634 ncoeff - 1 + k * 61] * mfactk;
3639 for (j = *ncoeff - 1; j >= i__1; --j) {
3640 bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
3642 for (nd = 1; nd <= i__2; ++nd) {
3643 pntcrb[nd] = pntcrb[nd] * *tparam + courbe[nd + j * courbe_dim1] *
3650 /* -------------------------------- The end -----------------------------
3658 //=======================================================================
3659 //function : AdvApp2Var_MathBase::mmeps1_
3661 //=======================================================================
3662 int AdvApp2Var_MathBase::mmeps1_(doublereal *epsilo)
3665 /* ***********************************************************************
3670 /* Extraction of EPS1 from COMMON MPRCSN. EPS1 is spatial zero */
3671 /* equal to 1.D-9 */
3675 /* MPRCSN,PRECISON,EPS1. */
3677 /* INPUT ARGUMENTS : */
3678 /* ------------------ */
3681 /* OUTPUT ARGUMENTS : */
3682 /* ------------------- */
3683 /* EPSILO : Value of EPS1 (spatial zero (10**-9)) */
3685 /* COMMONS USED : */
3686 /* ---------------- */
3688 /* REFERENCES CALLED : */
3689 /* ----------------------- */
3691 /* DESCRIPTION/NOTES/LIMITATIONS : */
3692 /* ----------------------------------- */
3693 /* EPS1 is ABSOLUTE spatial zero, so it is necessary */
3694 /* to use it whenever it is necessary to test if a variable */
3695 /* is null. For example, if the norm of a vector is lower than */
3696 /* EPS1, this vector is NULL ! (when one works in */
3697 /* REAL*8) It is absolutely not advised to test arguments */
3698 /* compared to EPS1**2. Taking into account the rounding errors inevitable */
3699 /* during calculations, this causes testing compared to 0.D0. */
3701 /* ***********************************************************************
3706 /* ***********************************************************************
3711 /* Gives tolerances of invalidity in stream */
3712 /* as well as limits of iterative processes */
3714 /* general context, modifiable by the user */
3718 /* PARAMETER , TOLERANCE */
3720 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3721 /* ----------------------------------- */
3722 /* INITIALISATION : profile , **VIA MPRFTX** at input in stream
3723 /* loading of default values of the profile in MPRFTX at input */
3724 /* in stream. They are preserved in local variables of MPRFTX */
3726 /* Reset of default values : MDFINT */
3727 /* Interactive modification by the user : MDBINT */
3729 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
3730 /* MEPSPB ... EPS3,EPS4 */
3731 /* MEPSLN ... EPS2, NITERM , NITERR */
3732 /* MEPSNR ... EPS2 , NITERM */
3733 /* MITERR ... NITERR */
3735 /* ***********************************************************************
3738 /* NITERM : max nb of iterations */
3739 /* NITERR : nb of rapid iterations */
3740 /* EPS1 : tolerance of 3D null distance */
3741 /* EPS2 : tolerance of parametric null distance */
3742 /* EPS3 : tolerance to avoid division by 0.. */
3743 /* EPS4 : angular tolerance */
3747 /* ***********************************************************************
3749 *epsilo = mmprcsn_.eps1;
3754 //=======================================================================
3755 //function : mmexthi_
3757 //=======================================================================
3758 int mmexthi_(integer *ndegre,
3762 /* System generated locals */
3765 /* Local variables */
3766 static integer iadd, ideb, ndeg2, nmod2, ii, ibb;
3769 /* **********************************************************************
3774 /* Extract of common LDGRTL the weight of formulas of */
3775 /* Gauss quadrature on all roots of Legendre polynoms of degree */
3776 /* NDEGRE defined on [-1,1]. */
3780 /* ALL, AB_SPECIFI::COMMON&, EXTRACTION, &WEIGHT, &GAUSS. */
3782 /* INPUT ARGUMENTS : */
3783 /* ------------------ */
3784 /* NDEGRE : Mathematic degree of Legendre polynom. It should have */
3785 /* 2 <= NDEGRE <= 61. */
3787 /* OUTPUT ARGUMENTS : */
3788 /* ------------------- */
3789 /* HWGAUS : The table of weights of Gauss quadrature formulas */
3790 /* relative to NDEGRE roots of a polynome de Legendre de */
3793 /* COMMONS UTILISES : */
3794 /* ---------------- */
3797 /* REFERENCES CALLED : */
3798 /* ----------------------- */
3800 /* DESCRIPTION/NOTES/LIMITATIONS : */
3801 /* ----------------------------------- */
3802 /* ATTENTION: The condition on NDEGRE ( 2 <= NDEGRE <= 61) is not */
3803 /* tested. The caller should make the test.
3805 /* Name of the routine */
3808 /* Common MLGDRTL: */
3809 /* This common includes POSITIVE roots of Legendre polynims */
3810 /* AND weights of Gauss quadrature formulas on all */
3811 /* POSITIVE roots of Legendre polynoms. */
3815 /* ***********************************************************************
3820 /* The common of Legendre roots. */
3826 /* DESCRIPTION/NOTES/LIMITATIONS : */
3827 /* ----------------------------------- */
3829 /* ***********************************************************************
3835 /* ROOTAB : Table of all roots of Legendre polynoms */
3836 /* within the interval [0,1]. They are ranked for the degrees increasing from */
3838 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
3839 /* The adressing is the same. */
3840 /* HI0TAB : Table of Legendre interpolators for root x=0 */
3841 /* of polynoms of UNEVEN degree. */
3842 /* RTLTB0 : Table of Li(uk) where uk are the roots of */
3843 /* Legendre polynom of EVEN degree. */
3844 /* RTLTB1 : Table of Li(uk) where uk are the roots of */
3845 /* Legendre polynom of UNEVEN degree. */
3848 /************************************************************************
3850 /* Parameter adjustments */
3854 ibb = AdvApp2Var_SysBase::mnfndeb_();
3856 AdvApp2Var_SysBase::mgenmsg_("MMEXTHI", 7L);
3859 ndeg2 = *ndegre / 2;
3860 nmod2 = *ndegre % 2;
3862 /* Address of Gauss weight associated to the 1st strictly */
3863 /* positive root of Legendre polynom of degree NDEGRE in MLGDRTL. */
3865 iadd = ndeg2 * (ndeg2 - 1) / 2 + 1;
3867 /* Index of the 1st HWGAUS element associated to the 1st strictly */
3868 /* positive root of Legendre polynom of degree NDEGRE. */
3870 ideb = (*ndegre + 1) / 2 + 1;
3872 /* Reading of weights associated to strictly positive roots. */
3875 for (ii = ideb; ii <= i__1; ++ii) {
3876 kpt = iadd + ii - ideb;
3877 hwgaus[ii] = mlgdrtl_.hiltab[kpt + nmod2 * 465 - 1];
3881 /* For strictly negative roots, the weight is the same. */
3882 /* i.e HW(1) = HW(NDEGRE), HW(2) = HW(NDEGRE-1), etc... */
3885 for (ii = 1; ii <= i__1; ++ii) {
3886 hwgaus[ii] = hwgaus[*ndegre + 1 - ii];
3890 /* Case of uneven NDEGRE, 0 is root of Legendre polynom, */
3891 /* associated Gauss weights are loaded. */
3894 hwgaus[ndeg2 + 1] = mlgdrtl_.hi0tab[ndeg2];
3897 /* --------------------------- The end ----------------------------------
3901 AdvApp2Var_SysBase::mgsomsg_("MMEXTHI", 7L);
3906 //=======================================================================
3907 //function : mmextrl_
3909 //=======================================================================
3910 int mmextrl_(integer *ndegre,
3913 /* System generated locals */
3916 /* Local variables */
3917 static integer iadd, ideb, ndeg2, nmod2, ii, ibb;
3921 /* **********************************************************************
3926 /* Extract of the Common LDGRTL of Legendre polynom roots */
3927 /* of degree NDEGRE defined on [-1,1]. */
3931 /* ALL, AB_SPECIFI::COMMON&, EXTRACTION, &ROOT, &LEGENDRE. */
3933 /* INPUT ARGUMENTS : */
3934 /* ------------------ */
3935 /* NDEGRE : Mathematic degree of Legendre polynom. */
3936 /* It is required to have 2 <= NDEGRE <= 61. */
3938 /* OUTPUT ARGUMENTS : */
3939 /* ------------------- */
3940 /* ROOTLG : The table of roots of Legendre polynom of degree */
3941 /* NDEGRE defined on [-1,1]. */
3943 /* COMMONS USED : */
3944 /* ---------------- */
3947 /* REFERENCES CALLED : */
3948 /* ----------------------- */
3950 /* DESCRIPTION/NOTES/LIMITATIONS : */
3951 /* ----------------------------------- */
3952 /* ATTENTION: Condition of NDEGRE ( 2 <= NDEGRE <= 61) is not */
3953 /* tested. The caller should make the test. */
3955 /* **********************************************************************
3959 /* Name of the routine */
3962 /* Common MLGDRTL: */
3963 /* This common includes POSITIVE roots of Legendre polynoms */
3964 /* AND the weight of Gauss quadrature formulas on all */
3965 /* POSITIVE roots of Legendre polynoms. */
3967 /* ***********************************************************************
3972 /* The common of Legendre roots. */
3979 /* ***********************************************************************
3982 /* ROOTAB : Table of all roots of Legendre polynoms */
3983 /* within the interval [0,1]. They are ranked for the degrees increasing from */
3985 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
3986 /* The adressing is the same. */
3987 /* HI0TAB : Table of Legendre interpolators for root x=0 */
3988 /* of polynoms of UNEVEN degree. */
3989 /* RTLTB0 : Table of Li(uk) where uk are the roots of */
3990 /* Legendre polynom of EVEN degree. */
3991 /* RTLTB1 : Table of Li(uk) where uk are the roots of */
3992 /* Legendre polynom of UNEVEN degree. */
3995 /************************************************************************
3997 /* Parameter adjustments */
4001 ibb = AdvApp2Var_SysBase::mnfndeb_();
4003 AdvApp2Var_SysBase::mgenmsg_("MMEXTRL", 7L);
4006 ndeg2 = *ndegre / 2;
4007 nmod2 = *ndegre % 2;
4009 /* Address of the 1st strictly positive root of Legendre polynom */
4010 /* of degree NDEGRE in MLGDRTL. */
4012 iadd = ndeg2 * (ndeg2 - 1) / 2 + 1;
4014 /* Indice, in ROOTLG, of the 1st strictly positive root */
4015 /* of Legendre polynom of degree NDEGRE. */
4017 ideb = (*ndegre + 1) / 2 + 1;
4019 /* Reading of strictly positive roots. */
4022 for (ii = ideb; ii <= i__1; ++ii) {
4023 kpt = iadd + ii - ideb;
4024 rootlg[ii] = mlgdrtl_.rootab[kpt + nmod2 * 465 - 1];
4028 /* Strictly negative roots are equal to positive roots
4030 /* to the sign i.e RT(1) = -RT(NDEGRE), RT(2) = -RT(NDEGRE-1), etc...
4034 for (ii = 1; ii <= i__1; ++ii) {
4035 rootlg[ii] = -rootlg[*ndegre + 1 - ii];
4039 /* Case NDEGRE uneven, 0 is root of Legendre polynom. */
4042 rootlg[ndeg2 + 1] = 0.;
4045 /* -------------------------------- THE END -----------------------------
4049 AdvApp2Var_SysBase::mgenmsg_("MMEXTRL", 7L);
4054 //=======================================================================
4055 //function : AdvApp2Var_MathBase::mmfmca8_
4057 //=======================================================================
4058 int AdvApp2Var_MathBase::mmfmca8_(integer *ndimen,
4068 /* System generated locals */
4069 integer tabini_dim1, tabini_dim2, tabini_offset, tabres_dim1, tabres_dim2,
4072 /* Local variables */
4073 static integer i__, j, k, ilong;
4077 /* **********************************************************************
4082 /* Expansion of a table containing only most important things into a */
4083 /* greater data table. */
4087 /* ALL, MATH_ACCES:: CARREAU&, DECOMPRESSION, &CARREAU */
4089 /* INPUT ARGUMENTS : */
4090 /* ------------------ */
4091 /* NDIMEN: Dimension of the workspace. */
4092 /* NCOEFU: Degree +1 of the table by u. */
4093 /* NCOEFV: Degree +1 of the table by v. */
4094 /* NDIMAX: Max dimension of the space. */
4095 /* NCFUMX: Max Degree +1 of the table by u. */
4096 /* NCFVMX: Max Degree +1 of the table by v. */
4097 /* TABINI: The table to be decompressed. */
4099 /* OUTPUT ARGUMENTS : */
4100 /* ------------------- */
4101 /* TABRES: Decompressed table. */
4103 /* COMMONS USED : */
4104 /* ---------------- */
4106 /* REFERENCES CALLED : */
4107 /* ----------------------- */
4109 /* DESCRIPTION/NOTES/LIMITATIONS : */
4110 /* ----------------------------------- */
4111 /* The following call : */
4113 /* CALL MMFMCA8(NDIMEN,NCOEFU,NCOEFV,NDIMAX,NCFUMX,NCFVMX,TABINI,TABINI)
4116 /* where TABINI is input/output argument, is possible provided */
4117 /* that the caller has declared TABINI in (NDIMAX,NCFUMX,NCFVMX) */
4119 /* ATTENTION : it is not checked that NDIMAX >= NDIMEN, */
4120 /* NCOEFU >= NCFMXU and NCOEFV >= NCFMXV. */
4122 /* **********************************************************************
4126 /* Parameter adjustments */
4127 tabini_dim1 = *ndimen;
4128 tabini_dim2 = *ncoefu;
4129 tabini_offset = tabini_dim1 * (tabini_dim2 + 1) + 1;
4130 tabini -= tabini_offset;
4131 tabres_dim1 = *ndimax;
4132 tabres_dim2 = *ncfumx;
4133 tabres_offset = tabres_dim1 * (tabres_dim2 + 1) + 1;
4134 tabres -= tabres_offset;
4137 if (*ndimax == *ndimen) {
4141 /* ----------------------- decompression NDIMAX<>NDIMEN -----------------
4144 for (k = *ncoefv; k >= 1; --k) {
4145 for (j = *ncoefu; j >= 1; --j) {
4146 for (i__ = *ndimen; i__ >= 1; --i__) {
4147 tabres[i__ + (j + k * tabres_dim2) * tabres_dim1] = tabini[
4148 i__ + (j + k * tabini_dim2) * tabini_dim1];
4157 /* ----------------------- decompression NDIMAX=NDIMEN ------------------
4161 if (*ncoefu == *ncfumx) {
4164 ilong = (*ndimen << 3) * *ncoefu;
4165 for (k = *ncoefv; k >= 1; --k) {
4166 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
4167 (char *)&tabini[(k * tabini_dim2 + 1) * tabini_dim1 + 1],
4168 (char *)&tabres[(k * tabres_dim2 + 1) * tabres_dim1 + 1]);
4173 /* ----------------- decompression NDIMAX=NDIMEN,NCOEFU=NCFUMX ----------
4177 ilong = (*ndimen << 3) * *ncoefu * *ncoefv;
4178 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
4179 (char *)&tabini[tabini_offset],
4180 (char *)&tabres[tabres_offset]);
4183 /* ---------------------------- The end ---------------------------------
4190 //=======================================================================
4191 //function : AdvApp2Var_MathBase::mmfmca9_
4193 //=======================================================================
4194 int AdvApp2Var_MathBase::mmfmca9_(integer *ndimax,
4204 /* System generated locals */
4205 integer tabini_dim1, tabini_dim2, tabini_offset, tabres_dim1, tabres_dim2,
4206 tabres_offset, i__1, i__2, i__3;
4208 /* Local variables */
4209 static integer i__, j, k, ilong;
4213 /* **********************************************************************
4218 /* Compression of a data table in a table */
4219 /* containing only the main data (the input table is not removed). */
4223 /* ALL, MATH_ACCES:: CARREAU&, COMPRESSION, &CARREAU */
4225 /* INPUT ARGUMENTS : */
4226 /* ------------------ */
4227 /* NDIMAX: Max dimension of the space. */
4228 /* NCFUMX: Max degree +1 of the table by u. */
4229 /* NCFVMX: Max degree +1 of the table by v. */
4230 /* NDIMEN: Dimension of the workspace. */
4231 /* NCOEFU: Degree +1 of the table by u. */
4232 /* NCOEFV: Degree +1 of the table by v. */
4233 /* TABINI: The table to compress. */
4235 /* OUTPUT ARGUMENTS : */
4236 /* ------------------- */
4237 /* TABRES: The compressed table. */
4239 /* COMMONS USED : */
4240 /* ---------------- */
4242 /* REFERENCES CALLED : */
4243 /* ----------------------- */
4245 /* DESCRIPTION/NOTES/LIMITATIONS : */
4246 /* ----------------------------------- */
4247 /* The following call : */
4249 /* CALL MMFMCA9(NDIMAX,NCFUMX,NCFVMX,NDIMEN,NCOEFU,NCOEFV,TABINI,TABINI)
4252 /* where TABINI is input/output argument, is possible provided */
4253 /* that the caller has checked that : */
4255 /* NDIMAX > NDIMEN, */
4256 /* or NDIMAX = NDIMEN and NCFUMX > NCOEFU */
4257 /* or NDIMAX = NDIMEN, NCFUMX = NCOEFU and NCFVMX > NCOEFV */
4259 /* These conditions are not tested in the program. */
4262 /* **********************************************************************
4266 /* Parameter adjustments */
4267 tabini_dim1 = *ndimax;
4268 tabini_dim2 = *ncfumx;
4269 tabini_offset = tabini_dim1 * (tabini_dim2 + 1) + 1;
4270 tabini -= tabini_offset;
4271 tabres_dim1 = *ndimen;
4272 tabres_dim2 = *ncoefu;
4273 tabres_offset = tabres_dim1 * (tabres_dim2 + 1) + 1;
4274 tabres -= tabres_offset;
4277 if (*ndimen == *ndimax) {
4281 /* ----------------------- Compression NDIMEN<>NDIMAX -------------------
4285 for (k = 1; k <= i__1; ++k) {
4287 for (j = 1; j <= i__2; ++j) {
4289 for (i__ = 1; i__ <= i__3; ++i__) {
4290 tabres[i__ + (j + k * tabres_dim2) * tabres_dim1] = tabini[
4291 i__ + (j + k * tabini_dim2) * tabini_dim1];
4300 /* ----------------------- Compression NDIMEN=NDIMAX --------------------
4304 if (*ncoefu == *ncfumx) {
4307 ilong = (*ndimen << 3) * *ncoefu;
4309 for (k = 1; k <= i__1; ++k) {
4310 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
4311 (char *)&tabini[(k * tabini_dim2 + 1) * tabini_dim1 + 1],
4312 (char *)&tabres[(k * tabres_dim2 + 1) * tabres_dim1 + 1]);
4317 /* ----------------- Compression NDIMEN=NDIMAX,NCOEFU=NCFUMX ------------
4321 ilong = (*ndimen << 3) * *ncoefu * *ncoefv;
4322 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
4323 (char *)&tabini[tabini_offset],
4324 (char *)&tabres[tabres_offset]);
4327 /* ---------------------------- The end ---------------------------------
4334 //=======================================================================
4335 //function : AdvApp2Var_MathBase::mmfmcar_
4337 //=======================================================================
4338 int AdvApp2Var_MathBase::mmfmcar_(integer *ndimen,
4351 static integer c__8 = 8;
4352 /* System generated locals */
4353 integer patold_dim1, patold_dim2, patnew_dim1, patnew_dim2,
4354 i__1, patold_offset,patnew_offset;
4356 /* Local variables */
4357 static doublereal tbaux[1];
4358 static integer ksize, numax, kk;
4359 static long int iofst;
4360 static integer ibb, ier;
4362 /* ***********************************************************************
4367 /* LIMITATION OF A SQUARE DEFINED ON (0,1)*(0,1) BETWEEN ISOS */
4368 /* UPARA1 AND UPARA2 (BY U) AND VPARA1 AND VPARA2 BY V. */
4372 /* LIMITATION , SQUARE , PARAMETER */
4374 /* INPUT ARGUMENTS : */
4375 /* ------------------ */
4376 /* NCOFMX: MAX NUMBER OF COEFF OF THE SQUARE BY U */
4377 /* NCOEFU: NUMBER OF COEFF OF THE SQUARE BY U */
4378 /* NCOEFV: NUMBER OF COEFF OF THE SQUARE BY V */
4379 /* PATOLD : THE SQUARE IS LIMITED BY UPARA1,UPARA2 AND VPARA1,VPARA2
4381 /* UPARA1 : LOWER LIMIT OF U */
4382 /* UPARA2 : UPPER LIMIT OF U */
4383 /* VPARA1 : LOWER LIMIT OF V */
4384 /* VPARA2 : UPPER LIMIT OF V */
4386 /* OUTPUT ARGUMENTS : */
4387 /* ------------------- */
4388 /* PATNEW : RELIMITED SQUARE, DEFINED ON (0,1)**2 */
4389 /* IERCOD : =10 COEFF NB TOO GREAT OR NULL */
4390 /* =13 PB IN THE DYNAMIC ALLOCATION */
4393 /* COMMONS USED : */
4394 /* ---------------- */
4396 /* DESCRIPTION/NOTES/LIMITATIONS : */
4397 /* ----------------------------------- */
4398 /* ---> The following call : */
4399 /* CALL MMFMCAR(NCOFMX,NCOEFU,NCOEFV,PATOLD,UPARA1,UPARA2,VPARA1,VPARA2
4402 /* where PATOLD is input/output argument is absolutely legal. */
4404 /* ---> The max number of coeff by u and v of PATOLD is 61 */
4406 /* ---> If NCOEFU < NCOFMX, the data is compressed by MMFMCA9 before
4407 /* limitation by v to get time during the execution */
4408 /* of MMARC41 that follows (the square is processed as a curve of
4410 /* dimension NDIMEN*NCOEFU possessing NCOEFV coefficients). */
4412 /* ***********************************************************************
4415 /* Name of the routine */
4418 /* Parameter adjustments */
4419 patnew_dim1 = *ndimen;
4420 patnew_dim2 = *ncofmx;
4421 patnew_offset = patnew_dim1 * (patnew_dim2 + 1) + 1;
4422 patnew -= patnew_offset;
4423 patold_dim1 = *ndimen;
4424 patold_dim2 = *ncofmx;
4425 patold_offset = patold_dim1 * (patold_dim2 + 1) + 1;
4426 patold -= patold_offset;
4429 ibb = AdvApp2Var_SysBase::mnfndeb_();
4431 AdvApp2Var_SysBase::mgenmsg_("MMFMCAR", 7L);
4436 /* **********************************************************************
4438 /* TEST OF COEFFICIENT NUMBERS */
4439 /* **********************************************************************
4442 if (*ncofmx < *ncoefu) {
4446 if (*ncoefu < 1 || *ncoefu > 61 || *ncoefv < 1 || *ncoefv > 61) {
4451 /* **********************************************************************
4453 /* CASE WHEN UPARA1=VPARA1=0 AND UPARA2=VPARA2=1 */
4454 /* **********************************************************************
4457 if (*upara1 == 0. && *upara2 == 1. && *vpara1 == 0. && *vpara2 == 1.) {
4458 ksize = (*ndimen << 3) * *ncofmx * *ncoefv;
4459 AdvApp2Var_SysBase::mcrfill_((integer *)&ksize,
4460 (char *)&patold[patold_offset],
4461 (char *)&patnew[patnew_offset]);
4465 /* **********************************************************************
4467 /* LIMITATION BY U */
4468 /* **********************************************************************
4471 if (*upara1 == 0. && *upara2 == 1.) {
4475 for (kk = 1; kk <= i__1; ++kk) {
4476 mmarc41_(ndimen, ndimen, ncoefu, &patold[(kk * patold_dim2 + 1) *
4477 patold_dim1 + 1], upara1, upara2, &patnew[(kk * patnew_dim2 +
4478 1) * patnew_dim1 + 1], iercod);
4482 /* **********************************************************************
4484 /* LIMITATION BY V */
4485 /* **********************************************************************
4489 if (*vpara1 == 0. && *vpara2 == 1.) {
4493 /* ----------- LIMITATION BY V (WITH COMPRESSION I.E. NCOEFU<NCOFMX) ----
4496 numax = *ndimen * *ncoefu;
4497 if (*ncofmx != *ncoefu) {
4498 /* ------------------------- Dynamic allocation -------------------
4500 ksize = *ndimen * *ncoefu * *ncoefv;
4501 AdvApp2Var_SysBase::mcrrqst_(&c__8, &ksize, tbaux, &iofst, &ier);
4506 /* --------------- Compression by (NDIMEN,NCOEFU,NCOEFV) ------------
4508 if (*upara1 == 0. && *upara2 == 1.) {
4509 AdvApp2Var_MathBase::mmfmca9_(ndimen,
4515 &patold[patold_offset],
4518 AdvApp2Var_MathBase::mmfmca9_(ndimen,
4524 &patnew[patnew_offset],
4527 /* ------------------------- Limitation by v ------------------------
4529 mmarc41_(&numax, &numax, ncoefv, &tbaux[iofst], vpara1, vpara2, &
4530 tbaux[iofst], iercod);
4531 /* --------------------- Expansion of TBAUX into PATNEW -------------
4533 AdvApp2Var_MathBase::mmfmca8_(ndimen, ncoefu, ncoefv, ndimen, ncofmx, ncoefv, &tbaux[iofst]
4534 , &patnew[patnew_offset]);
4537 /* -------- LIMITATION BY V (WITHOUT COMPRESSION I.E. NCOEFU=NCOFMX) ---
4541 if (*upara1 == 0. && *upara2 == 1.) {
4542 mmarc41_(&numax, &numax, ncoefv, &patold[patold_offset], vpara1,
4543 vpara2, &patnew[patnew_offset], iercod);
4545 mmarc41_(&numax, &numax, ncoefv, &patnew[patnew_offset], vpara1,
4546 vpara2, &patnew[patnew_offset], iercod);
4551 /* **********************************************************************
4554 /* **********************************************************************
4559 AdvApp2Var_SysBase::mcrdelt_(&c__8, &ksize, tbaux, &iofst, &ier);
4565 /* ------------------------------ The end -------------------------------
4570 AdvApp2Var_SysBase::maermsg_("MMFMCAR", iercod, 7L);
4573 AdvApp2Var_SysBase::mgsomsg_("MMFMCAR", 7L);
4579 //=======================================================================
4580 //function : AdvApp2Var_MathBase::mmfmcb5_
4582 //=======================================================================
4583 int AdvApp2Var_MathBase::mmfmcb5_(integer *isenmsc,
4594 /* System generated locals */
4595 integer courb1_dim1, courb1_offset, courb2_dim1, courb2_offset, i__1,
4598 /* Local variables */
4599 static integer i__, nboct, nd;
4602 /* **********************************************************************
4607 /* Reformating (and eventual compression/decompression) of curve */
4608 /* (ndim,.) by (.,ndim) and vice versa. */
4612 /* ALL , MATH_ACCES :: */
4613 /* COURBE&, REORGANISATION,COMPRESSION,INVERSION , &COURBE */
4615 /* INPUT ARGUMENTS : */
4616 /* -------------------- */
4617 /* ISENMSC : required direction of the transfer : */
4618 /* 1 : passage of (NDIMEN,.) ---> (.,NDIMEN) direction to AB
4620 /* -1 : passage of (.,NDIMEN) ---> (NDIMEN,.) direction to TS,T
4622 /* NDIMAX : format / dimension */
4623 /* NCF1MX : format by t of COURB1 */
4624 /* if ISENMSC= 1 : COURB1: The curve to be processed (NDIMAX,.) */
4625 /* NCOEFF : number of coeff of the curve */
4626 /* NCF2MX : format by t of COURB2 */
4627 /* NDIMEN : dimension of the curve and format of COURB2 */
4628 /* if ISENMSC=-1 : COURB2: The curve to be processed (.,NDIMEN) */
4630 /* OUTPUT ARGUMENTS : */
4631 /* --------------------- */
4632 /* if ISENMSC= 1 : COURB2: The resulting curve (.,NDIMEN) */
4633 /* if ISENMSC=-1 : COURB1: The resulting curve (NDIMAX,.) */
4635 /* COMMONS USED : */
4636 /* ------------------ */
4638 /* REFERENCES CALLED : */
4639 /* --------------------- */
4641 /* DESCRIPTION/NOTES/LIMITATIONS : */
4642 /* ----------------------------------- */
4643 /* allow to process the usual transfers as follows : */
4644 /* | ---- ISENMSC = 1 ---- | | ---- ISENMSC =-1 ----- | */
4645 /* TS (3,21) --> (21,3) AB ; AB (21,3) --> (3,21) TS */
4646 /* TS (3,21) --> (NU,3) AB ; AB (NU,3) --> (3,21) TS */
4647 /* (3,NU) --> (21,3) AB ; AB (21,3) --> (3,NU) */
4648 /* (3,NU) --> (NU,3) AB ; AB (NU,3) --> (3,NU) */
4650 /* ***********************************************************************
4654 /* Parameter adjustments */
4655 courb1_dim1 = *ndimax;
4656 courb1_offset = courb1_dim1 + 1;
4657 courb1 -= courb1_offset;
4658 courb2_dim1 = *ncf2mx;
4659 courb2_offset = courb2_dim1 + 1;
4660 courb2 -= courb2_offset;
4663 if (*ndimen > *ndimax || *ncoeff > *ncf1mx || *ncoeff > *ncf2mx) {
4667 if (*ndimen == 1 && *ncf1mx == *ncf2mx) {
4668 nboct = *ncf2mx << 3;
4669 if (*isenmsc == 1) {
4670 AdvApp2Var_SysBase::mcrfill_((integer *)&nboct,
4671 (char *)&courb1[courb1_offset],
4672 (char *)&courb2[courb2_offset]);
4674 if (*isenmsc == -1) {
4675 AdvApp2Var_SysBase::mcrfill_((integer *)&nboct,
4676 (char *)&courb2[courb2_offset],
4677 (char *)&courb1[courb1_offset]);
4684 if (*isenmsc == 1) {
4686 for (nd = 1; nd <= i__1; ++nd) {
4688 for (i__ = 1; i__ <= i__2; ++i__) {
4689 courb2[i__ + nd * courb2_dim1] = courb1[nd + i__ *
4695 } else if (*isenmsc == -1) {
4697 for (nd = 1; nd <= i__1; ++nd) {
4699 for (i__ = 1; i__ <= i__2; ++i__) {
4700 courb1[nd + i__ * courb1_dim1] = courb2[i__ + nd *
4712 /* ***********************************************************************
4720 AdvApp2Var_SysBase::maermsg_("MMFMCB5", iercod, 7L);
4725 //=======================================================================
4726 //function : AdvApp2Var_MathBase::mmfmtb1_
4728 //=======================================================================
4729 int AdvApp2Var_MathBase::mmfmtb1_(integer *maxsz1,
4739 static integer c__8 = 8;
4741 /* System generated locals */
4742 integer table1_dim1, table1_offset, table2_dim1, table2_offset, i__1,
4745 /* Local variables */
4746 static doublereal work[1];
4747 static integer ilong, isize, ii, jj, ier;
4748 static long int iofst,iipt, jjpt;
4751 /************************************************************************
4756 /* Inversion of elements of a rectangular table (T1(i,j) */
4757 /* loaded in T2(j,i)) */
4761 /* ALL, MATH_ACCES :: TABLEAU&, INVERSION, &TABLEAU */
4763 /* INPUT ARGUMENTS : */
4764 /* ------------------ */
4765 /* MAXSZ1: Max Nb of elements by the 1st dimension of TABLE1. */
4766 /* TABLE1: Table of reals by two dimensions. */
4767 /* ISIZE1: Nb of useful elements of TABLE1 on the 1st dimension */
4768 /* JSIZE1: Nb of useful elements of TABLE1 on the 2nd dimension */
4769 /* MAXSZ2: Nb max of elements by the 1st dimension of TABLE2. */
4771 /* OUTPUT ARGUMENTS : */
4772 /* ------------------- */
4773 /* TABLE2: Table of reals by two dimensions, containing the transposition
4774 /* of the rectangular table TABLE1. */
4775 /* ISIZE2: Nb of useful elements of TABLE2 on the 1st dimension */
4776 /* JSIZE2: Nb of useful elements of TABLE2 on the 2nd dimension */
4777 /* IERCOD: Erroe coder. */
4779 /* = 1, error in the dimension of tables */
4780 /* ether MAXSZ1 < ISIZE1 (table TABLE1 too small). */
4781 /* or MAXSZ2 < JSIZE1 (table TABLE2 too small). */
4783 /* COMMONS USED : */
4784 /* ---------------- */
4786 /* REFERENCES CALLED : */
4787 /* ---------------------- */
4789 /* DESCRIPTION/NOTES/LIMITATIONS : */
4790 /* ----------------------------------- */
4791 /* It is possible to use TABLE1 as input and output table i.e. */
4793 /* CALL MMFMTB1(MAXSZ1,TABLE1,ISIZE1,JSIZE1,MAXSZ2,TABLE1 */
4794 /* ,ISIZE2,JSIZE2,IERCOD) */
4797 /* **********************************************************************
4801 /* Parameter adjustments */
4802 table1_dim1 = *maxsz1;
4803 table1_offset = table1_dim1 + 1;
4804 table1 -= table1_offset;
4805 table2_dim1 = *maxsz2;
4806 table2_offset = table2_dim1 + 1;
4807 table2 -= table2_offset;
4811 if (*isize1 > *maxsz1 || *jsize1 > *maxsz2) {
4816 isize = *maxsz2 * *isize1;
4817 AdvApp2Var_SysBase::mcrrqst_(&c__8, &isize, work, &iofst, &ier);
4822 /* DO NOT BE AFRAID OF CRUSHING. */
4825 for (ii = 1; ii <= i__1; ++ii) {
4826 iipt = (ii - 1) * *maxsz2 + iofst;
4828 for (jj = 1; jj <= i__2; ++jj) {
4829 jjpt = iipt + (jj - 1);
4830 work[jjpt] = table1[ii + jj * table1_dim1];
4836 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
4837 (char *)&work[iofst],
4838 (char *)&table2[table2_offset]);
4840 /* -------------- The number of elements of TABLE2 is returned ------------
4849 /* ------------------------------- THE END ------------------------------
4851 /* --> Invalid input. */
4855 /* --> Pb of allocation. */
4862 AdvApp2Var_SysBase::mcrdelt_(&c__8, &isize, work, &iofst, &ier);
4870 //=======================================================================
4871 //function : AdvApp2Var_MathBase::mmgaus1_
4873 //=======================================================================
4874 int AdvApp2Var_MathBase::mmgaus1_(integer *ndimf,
4891 /* System generated locals */
4894 /* Local variables */
4895 static integer ndeg;
4896 static doublereal h__[20];
4898 static doublereal t, u[20], x;
4899 static integer idimf;
4900 static doublereal c1x, c2x;
4901 /* **********************************************************************
4907 /* Calculate the integral of function BFUNX passed in parameter */
4908 /* between limits XD and XF . */
4909 /* The function should be calculated for any value */
4910 /* of the variable in the given interval.. */
4911 /* The method GAUSS-LEGENDRE is used.
4912 /* For explications refer to the book : */
4913 /* Complements de mathematiques a l'usage des Ingenieurs de */
4914 /* l'electrotechnique et des telecommunications. */
4915 /* Par Andre ANGOT - Collection technique et scientifique du CNET
4918 /* The degree of LEGENDRE polynoms used is passed in parameter.
4922 /* INTEGRATION,LEGENDRE,GAUSS */
4924 /* INPUT ARGUMENTS : */
4925 /* ------------------ */
4927 /* NDIMF : Dimension of the function */
4928 /* BFUNX : Function to integrate passed as argument */
4929 /* Should be declared as EXTERNAL in the call routine. */
4930 /* SUBROUTINE BFUNX(NDIMF,X,VAL,IER) */
4932 /* K : Parameter determining the degree of the LEGENDRE polynom that
4934 /* can take a value between 0 and 10. */
4935 /* The degree of the polynom is equal to 4 k, that is 4, 8,
4937 /* 12, 16, 20, 24, 28, 32, 36 and 40. */
4938 /* If K is not correct, the degree is set to 40 directly.
4940 /* XD : Lower limit of the interval of integration. */
4941 /* XF : Upper limit of the interval of integration. */
4942 /* SAUX1 : Auxiliary table */
4943 /* SAUX2 : Auxiliary table */
4945 /* OUTPUT ARGUMENTS : */
4946 /* ------------------- */
4948 /* SOMME : Value of the integral */
4949 /* NITER : Number of iterations to be carried out. */
4950 /* It is equal to the degree of the polynom. */
4952 /* IER : Error code : */
4953 /* < 0 ==> Attention - Warning */
4954 /* = 0 ==> Everything is OK */
4955 /* > 0 ==> Critical error - Apply special processing */
4956 /* ==> Error in the calculation of BFUNX (return code */
4957 /* of this routine */
4959 /* If error => SUM = 0 */
4961 /* COMMONS USED : */
4962 /* ----------------- */
4966 /* REFERENCES CALLED : */
4967 /* ---------------------- */
4970 /* @ BFUNX MVGAUS0 */
4972 /* DESCRIPTION/NOTES/LIMITATIONS : */
4973 /* --------------------------------- */
4975 /* See the explanations detailed in the listing */
4976 /* Use of the GAUSS method (orthogonal polynoms) */
4977 /* The symmetry of roots of these polynomes is used */
4978 /* Depending on K, the degree of the interpolated polynom grows.
4980 /* If you wish to calculate the integral with a given precision, */
4981 /* loop on k varying from 1 to 10 and test the difference of 2
4983 /* consecutive iterations. Stop the loop if this difference is less that
4984 /* an epsilon value set to 10E-6 for example. */
4985 /* If S1 and S2 are 2 successive iterations, test following this example :
4988 /* AF=DABS(S1-S2) */
4990 /* If AS < 1 test if FS < eps otherwise test if AF/AS < eps
4992 /* -- ----- ----- */
4994 /************************************************************************
4997 /************************************************************************
5002 /* ****** General Initialization */
5004 /* Parameter adjustments */
5010 AdvApp2Var_SysBase::mvriraz_((integer *)ndimf,
5014 /* ****** Loading of coefficients U and H ** */
5015 /* -------------------------------------------- */
5017 mvgaus0_(k, u, h__, &ndeg, iercod);
5022 /* ****** C1X => Medium interval point [XD,XF] */
5023 /* ****** C2X => 1/2 amplitude interval [XD,XF] */
5025 c1x = (*xf + *xd) * .5;
5026 c2x = (*xf - *xd) * .5;
5028 /* ---------------------------------------- */
5029 /* ****** Integration for degree NDEG ** */
5030 /* ---------------------------------------- */
5033 for (j = 1; j <= i__1; ++j) {
5037 (*bfunx)(ndimf, &x, &saux1[1], iercod);
5043 (*bfunx)(ndimf, &x, &saux2[1], iercod);
5049 for (idimf = 1; idimf <= i__2; ++idimf) {
5050 somme[idimf] += h__[j - 1] * (saux1[idimf] + saux2[idimf]);
5057 for (idimf = 1; idimf <= i__1; ++idimf) {
5058 somme[idimf] *= c2x;
5061 /* ****** End of sub-program ** */
5067 //=======================================================================
5068 //function : mmherm0_
5070 //=======================================================================
5071 int mmherm0_(doublereal *debfin,
5074 static integer c__576 = 576;
5075 static integer c__6 = 6;
5078 /* System generated locals */
5082 /* Local variables */
5083 static doublereal amat[36] /* was [6][6] */;
5084 static integer iord[2];
5085 static doublereal prod;
5086 static integer iord1, iord2;
5087 static doublereal miden[36] /* was [6][6] */;
5088 static integer ncmat;
5089 static doublereal epspi, d1, d2;
5090 static integer ii, jj, pp, ncf;
5091 static doublereal cof[6];
5092 static integer iof[2], ier;
5093 static doublereal mat[36] /* was [6][6] */;
5095 static doublereal abid[72] /* was [12][6] */;
5096 /* ***********************************************************************
5101 /* INIT OF COEFFS. OF POLYNOMS OF HERMIT INTERPOLATION */
5105 /* MATH_ACCES :: HERMITE */
5107 /* INPUT ARGUMENTS */
5108 /* -------------------- */
5109 /* DEBFIN : PARAMETERS DEFINING THE CONSTRAINTS */
5110 /* DEBFIN(1) : FIRST PARAMETER */
5111 /* DEBFIN(2) : SECOND PARAMETER */
5113 /* ONE SHOULD HAVE: */
5114 /* ABS (DEBFIN(I)) < 100 */
5116 /* (ABS(DEBFIN(1)+ABS(DEBFIN(2))) > 1/100 */
5117 /* (for overflows) */
5119 /* ABS(DEBFIN(2)-DEBFIN(1)) / (ABS(DEBFIN(1)+ABS(DEBFIN(2))) > 1/100
5121 /* (for the conditioning) */
5124 /* OUTPUT ARGUMENTS : */
5125 /* --------------------- */
5127 /* IERCOD : Error code : 0 : O.K. */
5128 /* 1 : value of DEBFIN */
5129 /* are unreasonable */
5130 /* -1 : init was already done */
5131 /* (OK but no processing) */
5133 /* COMMONS USED : */
5134 /* ------------------ */
5136 /* REFERENCES CALLED : */
5137 /* ---------------------- */
5140 /* DESCRIPTION/NOTES/LIMITATIONS : */
5141 /* ----------------------------------- */
5143 /* This program initializes the coefficients of Hermit polynoms */
5144 /* that are read later by MMHERM1 */
5145 /* ***********************************************************************
5150 /* **********************************************************************
5155 /* Used to STORE coefficients of Hermit interpolation polynoms
5161 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5162 /* ----------------------------------- */
5164 /* The coefficients of hermit polynoms are calculated by */
5165 /* the routine MMHERM0 and read by the routine MMHERM1 */
5167 /* **********************************************************************
5174 /* NBCOEF is the size of CMHERM (see below) */
5175 /* ***********************************************************************
5184 /* ***********************************************************************
5187 /* ***********************************************************************
5191 /* Parameter adjustments */
5195 d1 = abs(debfin[1]);
5196 if (d1 > (float)100.) {
5200 d2 = abs(debfin[2]);
5201 if (d2 > (float)100.) {
5206 if (d2 < (float).01) {
5210 d1 = (d__1 = debfin[2] - debfin[1], abs(d__1));
5211 if (d1 / d2 < (float).01) {
5216 /* ***********************************************************************
5218 /* Initialization */
5219 /* ***********************************************************************
5227 /* ***********************************************************************
5230 /* IS IT ALREADY INITIALIZED ? */
5232 d1 = abs(debfin[1]) + abs(debfin[2]);
5235 if (debfin[1] != mmcmher_.tdebut) {
5238 if (debfin[2] != mmcmher_.tfinal) {
5241 if (d1 != mmcmher_.verifi) {
5249 /* ***********************************************************************
5252 /* ***********************************************************************
5258 /* Init. matrix identity : */
5261 AdvApp2Var_SysBase::mvriraz_((integer *)&ncmat,
5264 for (ii = 1; ii <= 6; ++ii) {
5265 miden[ii + ii * 6 - 7] = 1.;
5271 /* Init to 0 of table CMHERM */
5273 AdvApp2Var_SysBase::mvriraz_((integer *)&c__576, (char *)mmcmher_.cmherm);
5275 /* Calculation by solution of linear systems */
5277 for (iord1 = -1; iord1 <= 2; ++iord1) {
5278 for (iord2 = -1; iord2 <= 2; ++iord2) {
5285 iof[1] = iord[0] + 1;
5288 ncf = iord[0] + iord[1] + 2;
5290 /* Calculate matrix MAT to invert: */
5292 for (cot = 1; cot <= 2; ++cot) {
5295 if (iord[cot - 1] > -1) {
5298 for (jj = 1; jj <= i__1; ++jj) {
5304 i__1 = iord[cot - 1] + 1;
5305 for (pp = 1; pp <= i__1; ++pp) {
5307 ii = pp + iof[cot - 1];
5312 for (jj = 1; jj <= i__2; ++jj) {
5313 mat[ii + jj * 6 - 7] = (float)0.;
5318 for (jj = pp; jj <= i__2; ++jj) {
5320 /* everything is done in these 3 lines
5323 mat[ii + jj * 6 - 7] = cof[jj - 1] * prod;
5324 cof[jj - 1] *= jj - pp;
5325 prod *= debfin[cot];
5338 AdvApp2Var_MathBase::mmmrslwd_(&c__6, &ncf, &ncf, mat, miden, &epspi, abid, amat, &
5345 for (cot = 1; cot <= 2; ++cot) {
5346 i__1 = iord[cot - 1] + 1;
5347 for (pp = 1; pp <= i__1; ++pp) {
5349 for (ii = 1; ii <= i__2; ++ii) {
5350 mmcmher_.cmherm[ii + (pp + (cot + ((iord1 + (iord2 <<
5351 2)) << 1)) * 3) * 6 + 155] = amat[ii + (pp +
5352 iof[cot - 1]) * 6 - 7];
5365 /* ***********************************************************************
5368 /* The initialized flag is located: */
5370 mmcmher_.tdebut = debfin[1];
5371 mmcmher_.tfinal = debfin[2];
5373 d1 = abs(debfin[1]) + abs(debfin[2]);
5374 mmcmher_.verifi = d1 * 16111959;
5377 /* ***********************************************************************
5382 /* ***********************************************************************
5393 /* ***********************************************************************
5398 AdvApp2Var_SysBase::maermsg_("MMHERM0", iercod, 7L);
5400 /* ***********************************************************************
5405 //=======================================================================
5406 //function : mmherm1_
5408 //=======================================================================
5409 int mmherm1_(doublereal *debfin,
5415 /* System generated locals */
5416 integer hermit_dim1, hermit_dim2, hermit_offset;
5418 /* Local variables */
5419 static integer nbval;
5420 static doublereal d1;
5423 /* ***********************************************************************
5428 /* reading of coeffs. of HERMIT interpolation polynoms */
5432 /* MATH_ACCES :: HERMIT */
5434 /* INPUT ARGUMENTS : */
5435 /* -------------------- */
5436 /* DEBFIN : PARAMETES DEFINING THE CONSTRAINTS */
5437 /* DEBFIN(1) : FIRST PARAMETER */
5438 /* DEBFIN(2) : SECOND PARAMETER */
5440 /* Should be equal to the corresponding arguments during the */
5441 /* last call to MMHERM0 for the initialization of coeffs. */
5443 /* ORDRMX : indicates the dimensioning of HERMIT: */
5444 /* there is no choice : ORDRMX should be equal to the value */
5445 /* of PARAMETER IORDMX of INCLUDE MMCMHER, or 2 for the moment */
5447 /* IORDRE (2) : Orders of constraints in each corresponding parameter DEBFIN(I)
5448 /* should be between -1 (no constraints) and ORDRMX. */
5451 /* OUTPUT ARGUMENTS : */
5452 /* --------------------- */
5454 /* HERMIT : HERMIT(1:IORDRE(1)+IORDRE(2)+2, j, cote) are the */
5455 /* coefficients in the canonic base of Hermit polynom */
5456 /* corresponding to orders IORDRE with parameters DEBFIN for */
5457 /* the constraint of order j on DEBFIN(cote). j is between 0 and IORDRE(cote). */
5460 /* IERCOD : Error code : */
5461 /* -1: O.K but necessary to reinitialize the coefficients */
5462 /* (info for optimization) */
5464 /* 1 : Error in MMHERM0 */
5465 /* 2 : arguments invalid */
5467 /* COMMONS USED : */
5468 /* ------------------ */
5470 /* REFERENCES CALLED : */
5471 /* ---------------------- */
5474 /* DESCRIPTION/NOTES/LIMITATIONS : */
5475 /* ----------------------------------- */
5477 /* This program reads coefficients of Hermit polynoms */
5478 /* that were earlier initialized by MMHERM0 */
5480 /* PMN : initialisation is no more done by the caller. */
5483 /* ***********************************************************************
5488 /* **********************************************************************
5493 /* Serves to STORE the coefficients of Hermit interpolation polynoms
5499 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5500 /* ----------------------------------- */
5502 /* the coefficients of Hetmit polynoms are calculated by */
5503 /* routine MMHERM0 and read by routine MMHERM1 */
5506 /* **********************************************************************
5513 /* NBCOEF is the size of CMHERM (see lower) */
5517 /* ***********************************************************************
5524 /* ***********************************************************************
5526 /* Initializations */
5527 /* ***********************************************************************
5530 /* Parameter adjustments */
5532 hermit_dim1 = (*ordrmx << 1) + 2;
5533 hermit_dim2 = *ordrmx + 1;
5534 hermit_offset = hermit_dim1 * hermit_dim2 + 1;
5535 hermit -= hermit_offset;
5542 /* ***********************************************************************
5545 /* ***********************************************************************
5553 for (cot = 1; cot <= 2; ++cot) {
5554 if (iordre[cot] < -1) {
5557 if (iordre[cot] > *ordrmx) {
5564 /* IS-IT CORRECTLY INITIALIZED ? */
5566 d1 = abs(debfin[1]) + abs(debfin[2]);
5569 /* OTHERWISE IT IS INITIALIZED */
5571 if (debfin[1] != mmcmher_.tdebut || debfin[2] != mmcmher_.tfinal || d1
5572 != mmcmher_.verifi) {
5574 mmherm0_(&debfin[1], iercod);
5581 /* ***********************************************************************
5584 /* ***********************************************************************
5589 AdvApp2Var_SysBase::msrfill_(&nbval, &mmcmher_.cmherm[((((iordre[1] + (iordre[2] << 2)) << 1)
5590 + 1) * 3 + 1) * 6 + 156], &hermit[hermit_offset]);
5592 /* ***********************************************************************
5597 /* ***********************************************************************
5608 /* ***********************************************************************
5613 AdvApp2Var_SysBase::maermsg_("MMHERM1", iercod, 7L);
5615 /* ***********************************************************************
5620 //=======================================================================
5621 //function : AdvApp2Var_MathBase::mmhjcan_
5623 //=======================================================================
5624 int AdvApp2Var_MathBase::mmhjcan_(integer *ndimen,
5635 static integer c__2 = 2;
5636 static integer c__21 = 21;
5637 /* System generated locals */
5638 integer tcbold_dim1, tcbold_dim2, tcbold_offset, tcbnew_dim1, tcbnew_dim2,
5639 tcbnew_offset, i__1, i__2, i__3, i__4, i__5;
5642 /* Local variables */
5643 static logical ldbg;
5644 static integer ndeg;
5645 static doublereal taux1[21];
5646 static integer d__, e, i__, k;
5647 static doublereal mfact;
5648 static integer ncoeff;
5649 static doublereal tjacap[21];
5650 static integer iordre[2];
5651 static doublereal hermit[36]/* was [6][3][2] */, ctenor, bornes[2];
5653 static integer aux1, aux2;
5655 /* ***********************************************************************
5660 /* CONVERSION OF TABLE TCBOLD OF POLYNOMIAL CURVE COEFFICIENTS */
5661 /* EXPRESSED IN HERMIT JACOBI BASE, INTO A */
5662 /* TABLE OF COEFFICIENTS TCBNEW OF COURVES EXPRESSED IN THE CANONIC BASE */
5666 /* CANNONIC, HERMIT, JACCOBI */
5668 /* INPUT ARGUMENTS : */
5669 /* -------------------- */
5670 /* ORDHER : ORDER OF HERMIT POLYNOMS OR ORDER OF CONTINUITY */
5671 /* NCOEFS : NUMBER OF COEFFICIENTS OF A POLYNOMIAL CURVE */
5672 /* FOR ONE OF ITS NDIM COMPONENTS;(DEGREE+1 OF THE CURVE)
5674 /* NDIM : DIMENSION OF THE CURVE */
5675 /* CBHEJA : TABLE OF COEFFICIENTS OF THE CURVE IN THE BASE */
5677 /* (H(0,-1),..,H(ORDHER,-1),H(0,1),..,H(ORDHER,1), */
5678 /* JA(ORDHER+1,2*ORDHER+2),....,JA(ORDHER+1,NCOEFS-1) */
5680 /* OUTPUT ARGUMENTS : */
5681 /* --------------------- */
5682 /* CBRCAN : TABLE OF COEFFICIENTS OF THE CURVE IN THE CANONIC BASE */
5685 /* COMMONS USED : */
5686 /* ------------------ */
5689 /* REFERENCES CALLED : */
5690 /* --------------------- */
5693 /* ***********************************************************************
5697 /* ***********************************************************************
5702 /* Providesinteger constants from 0 to 1000 */
5708 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5709 /* ----------------------------------- */
5711 /* ***********************************************************************
5715 /* ***********************************************************************
5721 /* ***********************************************************************
5723 /* INITIALIZATION */
5724 /* ***********************************************************************
5727 /* Parameter adjustments */
5729 tcbnew_dim1 = *ndimen;
5730 tcbnew_dim2 = *ncflim;
5731 tcbnew_offset = tcbnew_dim1 * (tcbnew_dim2 + 1) + 1;
5732 tcbnew -= tcbnew_offset;
5733 tcbold_dim1 = *ndimen;
5734 tcbold_dim2 = *ncflim;
5735 tcbold_offset = tcbold_dim1 * (tcbold_dim2 + 1) + 1;
5736 tcbold -= tcbold_offset;
5739 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
5741 AdvApp2Var_SysBase::mgenmsg_("MMHJCAN", 7L);
5748 /* ***********************************************************************
5751 /* ***********************************************************************
5761 /* CALCULATION OF HERMIT POLYNOMS IN THE CANONIC BASE ON (-1,1) */
5764 iordre[0] = *orcont;
5765 iordre[1] = *orcont;
5766 mmherm1_(bornes, &c__2, iordre, hermit, &ier);
5776 for (e = 1; e <= i__1; ++e) {
5778 ctenor = (tdecop[e] - tdecop[e - 1]) / 2;
5786 for (d__ = 1; d__ <= i__2; ++d__) {
5788 /* CONVERSION OF THE COEFFICIENTS OF THE PART OF THE CURVE EXPRESSED */
5789 /* IN HERMIT BASE, INTO THE CANONIC BASE */
5791 AdvApp2Var_SysBase::mvriraz_((integer *)&ncoeff, (char *)taux1);
5794 for (k = 1; k <= i__3; ++k) {
5796 for (i__ = 1; i__ <= i__4; ++i__) {
5798 mfact = AdvApp2Var_MathBase::pow__di(&ctenor, &i__5);
5799 taux1[k - 1] += (tcbold[d__ + (i__ + e * tcbold_dim2) *
5800 tcbold_dim1] * hermit[k + (i__ + 2) * 6 - 19] +
5801 tcbold[d__ + (i__ + aux1 + e * tcbold_dim2) *
5802 tcbold_dim1] * hermit[k + (i__ + 5) * 6 - 19]) *
5809 for (i__ = aux2 + 1; i__ <= i__3; ++i__) {
5810 taux1[i__ - 1] = tcbold[d__ + (i__ + e * tcbold_dim2) *
5814 /* CONVERSION OF THE COEFFICIENTS OF THE PART OF THE CURVE EXPRESSED */
5815 /* IN CANONIC-JACOBI BASE, INTO THE CANONIC BASE */
5819 AdvApp2Var_MathBase::mmapcmp_(&minombr_.nbr[1], &c__21, &ncoeff, taux1, tjacap);
5820 AdvApp2Var_MathBase::mmjacan_(orcont, &ndeg, tjacap, taux1);
5822 /* RECOPY THE COEFS RESULTING FROM THE CONVERSION IN THE TABLE */
5826 for (i__ = 1; i__ <= i__3; ++i__) {
5827 tcbnew[d__ + (i__ + e * tcbnew_dim2) * tcbnew_dim1] = taux1[
5836 /* ***********************************************************************
5838 /* PROCESSING OF ERRORS */
5839 /* ***********************************************************************
5849 /* ***********************************************************************
5851 /* RETURN CALLING PROGRAM */
5852 /* ***********************************************************************
5857 AdvApp2Var_SysBase::maermsg_("MMHJCAN", iercod, 7L);
5859 AdvApp2Var_SysBase::mgsomsg_("MMHJCAN", 7L);
5864 //=======================================================================
5865 //function : AdvApp2Var_MathBase::mminltt_
5867 //=======================================================================
5868 int AdvApp2Var_MathBase::mminltt_(integer *ncolmx,
5874 doublereal *,//epseg,
5877 /* System generated locals */
5878 integer tabtri_dim1, tabtri_offset, i__1, i__2;
5880 /* Local variables */
5881 static logical idbg;
5882 static integer icol, ilgn, nlgn, noct, inser;
5883 static doublereal epsega;
5886 /* ***********************************************************************
5891 /* . Insert a line in a table parsed without redundance */
5895 /* TOUS,MATH_ACCES :: TABLEAU&,INSERTION,&TABLEAU */
5897 /* INPUT ARGUMENTS : */
5898 /* -------------------- */
5899 /* . NCOLMX : Number of columns in the table */
5900 /* . NLGNMX : Number of lines in the table */
5901 /* . TABTRI : Table parsed by lines without redundances */
5902 /* . NBRCOL : Number of columns used */
5903 /* . NBRLGN : Number of lines used */
5904 /* . AJOUTE : Line to be added */
5905 /* . EPSEGA : Epsilon to test the redundance */
5907 /* OUTPUT ARGUMENTS : */
5908 /* --------------------- */
5909 /* . TABTRI : Table parsed by lines without redundances */
5910 /* . NBRLGN : Number of lines used */
5911 /* . IERCOD : 0 -> No problem */
5912 /* 1 -> The table is full */
5914 /* COMMONS USED : */
5915 /* ------------------ */
5917 /* REFERENCES CALLED : */
5918 /* --------------------- */
5920 /* DESCRIPTION/NOTES/LIMITATIONS : */
5921 /* ----------------------------------- */
5922 /* . The line is inserted only if there is no line with all
5924 /* elements equl to those which are planned to be insered, to epsilon. */
5926 /* . Level of de debug = 3 */
5930 /* DECLARATIONS , CONTROL OF INPUT ARGUMENTS , INITIALIZATION */
5931 /* ***********************************************************************
5934 /* --- Parameters */
5940 /* --- Local variables */
5945 /* Parameter adjustments */
5946 tabtri_dim1 = *ncolmx;
5947 tabtri_offset = tabtri_dim1 + 1;
5948 tabtri -= tabtri_offset;
5952 ibb = AdvApp2Var_SysBase::mnfndeb_();
5955 AdvApp2Var_SysBase::mgenmsg_("MMINLTT", 7L);
5958 /* --- Control arguments */
5960 if (*nbrlgn >= *nlgnmx) {
5964 /* -------------------- */
5965 /* *** INITIALIZATION */
5966 /* -------------------- */
5970 /* ---------------------------- */
5971 /* *** SEARCH OF REDUNDANCE */
5972 /* ---------------------------- */
5975 for (ilgn = 1; ilgn <= i__1; ++ilgn) {
5976 if (tabtri[ilgn * tabtri_dim1 + 1] >= ajoute[1] - epsega) {
5977 if (tabtri[ilgn * tabtri_dim1 + 1] <= ajoute[1] + epsega) {
5979 for (icol = 1; icol <= i__2; ++icol) {
5980 if (tabtri[icol + ilgn * tabtri_dim1] < ajoute[icol] -
5981 epsega || tabtri[icol + ilgn * tabtri_dim1] >
5982 ajoute[icol] + epsega) {
5996 /* ----------------------------------- */
5997 /* *** SEARCH OF THE INSERTION POINT */
5998 /* ----------------------------------- */
6003 for (ilgn = 1; ilgn <= i__1; ++ilgn) {
6005 for (icol = 1; icol <= i__2; ++icol) {
6006 if (tabtri[icol + ilgn * tabtri_dim1] < ajoute[icol]) {
6009 if (tabtri[icol + ilgn * tabtri_dim1] > ajoute[icol]) {
6020 /* -------------- */
6022 /* -------------- */
6029 /* --- Shift lower */
6031 nlgn = *nbrlgn - inser;
6033 noct = (*ncolmx << 3) * nlgn;
6034 AdvApp2Var_SysBase::mcrfill_((integer *)&noct,
6035 (char *)&tabtri[inser * tabtri_dim1 + 1],
6036 (char *)&tabtri[(inser + 1)* tabtri_dim1 + 1]);
6041 noct = *nbrcol << 3;
6042 AdvApp2Var_SysBase::mcrfill_((integer *)&noct,
6044 (char *)&tabtri[inser * tabtri_dim1 + 1]);
6048 /* ******************************************************************** */
6049 /* OUTPUT ERROR , RETURN CALLING PROGRAM , MESSAGES */
6050 /* ******************************************************************** */
6052 /* --- The table is already full */
6061 AdvApp2Var_SysBase::maermsg_("MMINLTT", iercod, 7L);
6064 AdvApp2Var_SysBase::mgsomsg_("MMINLTT", 7L);
6069 //=======================================================================
6070 //function : AdvApp2Var_MathBase::mmjacan_
6072 //=======================================================================
6073 int AdvApp2Var_MathBase::mmjacan_(integer *ideriv,
6078 /* System generated locals */
6079 integer poljac_dim1, i__1, i__2;
6081 /* Local variables */
6082 static integer iptt, i__, j, ibb;
6083 static doublereal bid;
6085 /* ***********************************************************************
6090 /* Routine of transfer of Jacobi normalized to canonic [-1,1], */
6091 /* the tables are ranked by even, then by uneven degree. */
6095 /* LEGENDRE,JACOBI,PASSAGE. */
6097 /* INPUT ARGUMENTS : */
6098 /* ------------------ */
6099 /* IDERIV : Order of Jacobi between -1 and 2. */
6100 /* NDEG : The true degree of the polynom. */
6101 /* POLJAC : The polynom in the Jacobi base. */
6103 /* OUTPUT ARGUMENTS : */
6104 /* ------------------- */
6105 /* POLCAN : The curve expressed in the canonic base [-1,1]. */
6107 /* COMMONS USED : */
6108 /* ---------------- */
6110 /* REFERENCES CALLED : */
6111 /* ----------------------- */
6113 /* DESCRIPTION/NOTES/LIMITATIONS : */
6114 /* ----------------------------------- */
6117 /* ***********************************************************************
6120 /* Name of the routine */
6122 /* Matrices of conversion */
6125 /* ***********************************************************************
6130 /* MATRIX OF TRANSFORMATION OF LEGENDRE BASE */
6136 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
6137 /* ----------------------------------- */
6140 /* ***********************************************************************
6145 /* Legendre common / Restricted Casteljau. */
6147 /* 0:1 0 Concerns the even terms, 1 the uneven terms. */
6148 /* CANPLG : Matrix of passage to canonic from Jacobi with calculated parities */
6149 /* PLGCAN : Matrix of passage from Jacobi to canonic with calculated parities */
6152 /* ***********************************************************************
6155 /* Parameter adjustments */
6156 poljac_dim1 = *ndeg / 2 + 1;
6159 ibb = AdvApp2Var_SysBase::mnfndeb_();
6161 AdvApp2Var_SysBase::mgenmsg_("MMJACAN", 7L);
6164 /* ----------------- Expression of terms of even degree ----------------
6168 for (i__ = 0; i__ <= i__1; ++i__) {
6170 iptt = i__ * 31 - (i__ + 1) * i__ / 2 + 1;
6172 for (j = i__; j <= i__2; ++j) {
6173 bid += mmjcobi_.plgcan[iptt + j + *ideriv * 992 + 991] * poljac[
6177 polcan[i__ * 2] = bid;
6181 /* --------------- Expression of terms of uneven degree ----------------
6188 i__1 = (*ndeg - 1) / 2;
6189 for (i__ = 0; i__ <= i__1; ++i__) {
6191 iptt = i__ * 31 - (i__ + 1) * i__ / 2 + 1;
6192 i__2 = (*ndeg - 1) / 2;
6193 for (j = i__; j <= i__2; ++j) {
6194 bid += mmjcobi_.plgcan[iptt + j + ((*ideriv << 1) + 1) * 496 +
6195 991] * poljac[j + poljac_dim1];
6198 polcan[(i__ << 1) + 1] = bid;
6202 /* -------------------------------- The end -----------------------------
6207 AdvApp2Var_SysBase::mgsomsg_("MMJACAN", 7L);
6212 //=======================================================================
6213 //function : AdvApp2Var_MathBase::mmjaccv_
6215 //=======================================================================
6216 int AdvApp2Var_MathBase::mmjaccv_(integer *ncoef,
6224 /* Initialized data */
6226 static char nomprg[8+1] = "MMJACCV ";
6228 /* System generated locals */
6229 integer crvlgd_dim1, crvlgd_offset, crvcan_dim1, crvcan_offset,
6230 polaux_dim1, i__1, i__2;
6232 /* Local variables */
6233 static integer ndeg, i__, nd, ii, ibb;
6235 /* ***********************************************************************
6240 /* Passage from the normalized Jacobi base to the canonic base. */
6244 /* SMOOTHING, BASE, LEGENDRE */
6247 /* INPUT ARGUMENTS : */
6248 /* ------------------ */
6249 /* NDIM: Space Dimension. */
6250 /* NCOEF: Degree +1 of the polynom. */
6251 /* IDER: Order of Jacobi polynoms. */
6252 /* CRVLGD : Curve in the base of Jacobi. */
6254 /* OUTPUT ARGUMENTS : */
6255 /* ------------------- */
6256 /* POLAUX : Auxilliary space. */
6257 /* CRVCAN : The curve in the canonic base [-1,1] */
6259 /* COMMONS USED : */
6260 /* ---------------- */
6262 /* REFERENCES CALLED : */
6263 /* ----------------------- */
6265 /* DESCRIPTION/NOTES/LIMITATIONS : */
6266 /* ----------------------------------- */
6269 /* *********************************************************************
6272 /* Name of the routine */
6273 /* Parameter adjustments */
6274 polaux_dim1 = (*ncoef - 1) / 2 + 1;
6275 crvcan_dim1 = *ncoef - 1 + 1;
6276 crvcan_offset = crvcan_dim1;
6277 crvcan -= crvcan_offset;
6278 crvlgd_dim1 = *ncoef - 1 + 1;
6279 crvlgd_offset = crvlgd_dim1;
6280 crvlgd -= crvlgd_offset;
6284 ibb = AdvApp2Var_SysBase::mnfndeb_();
6286 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
6292 for (nd = 1; nd <= i__1; ++nd) {
6293 /* Loading of the auxilliary table. */
6296 for (i__ = 0; i__ <= i__2; ++i__) {
6297 polaux[i__] = crvlgd[ii + nd * crvlgd_dim1];
6304 i__2 = (ndeg - 1) / 2;
6305 for (i__ = 0; i__ <= i__2; ++i__) {
6306 polaux[i__ + polaux_dim1] = crvlgd[ii + nd * crvlgd_dim1];
6311 /* Call the routine of base change. */
6312 AdvApp2Var_MathBase::mmjacan_(ider, &ndeg, polaux, &crvcan[nd * crvcan_dim1]);
6321 //=======================================================================
6322 //function : mmloncv_
6324 //=======================================================================
6325 int mmloncv_(integer *ndimax,
6335 /* Initialized data */
6337 static integer kgar = 0;
6339 /* System generated locals */
6340 integer courbe_dim1, courbe_offset, i__1, i__2;
6342 /* Local variables */
6343 static doublereal tran;
6344 static integer ngaus;
6345 static doublereal c1, c2, d1, d2, wgaus[20], uroot[20], x1, x2, dd;
6346 static integer ii, jj, kk;
6347 static doublereal som;
6348 static doublereal der1, der2;
6353 /* **********************************************************************
6356 /* FUNCTION : Length of an arc of curve on a given interval */
6357 /* ---------- for a function the mathematic representation */
6358 /* which of is a multidimensional polynom. */
6359 /* The polynom is a set of polynoms the coefficients which of are ranked
6360 /* in a table with 2 indices, each line relative to 1 polynom. */
6361 /* The polynom is defined by its coefficients ordered by increasing
6362 * power of the variable. */
6363 /* All polynoms have the same number of coefficients (and the same degree). */
6365 /* KEYWORDS : LENGTH, CURVE */
6368 /* INPUT ARGUMENTS : */
6369 /* -------------------- */
6371 /* NDIMAX : Max number of lines of tables (max number of polynoms). */
6372 /* NDIMEN : Dimension of the polynom (Nomber of polynoms). */
6373 /* NCOEFF : Number of coefficients of the polynom (no limitation) */
6374 /* This is degree + 1 */
6375 /* COURBE : Coefficients of the polynom ordered by increasing power */
6376 /* Dimension to (NDIMAX,NCOEFF). */
6377 /* TDEBUT : Lower limit of integration for length calculation. */
6378 /* TFINAL : Upper limit of integration for length calculation. */
6380 /* OUTPUT ARGUMENTS : */
6381 /* --------------------- */
6382 /* XLONGC : Length of arc of curve */
6384 /* IERCOD : Error code : */
6385 /* = 0 ==> All is OK */
6386 /* = 1 ==> NDIMEN or NCOEFF negative or null */
6387 /* = 2 ==> Pb loading Legendre roots and Gauss weight */
6390 /* If error => XLONGC = 0 */
6392 /* COMMONS USED : */
6393 /* ------------------ */
6397 /* REFERENCES CALLED : */
6398 /* ---------------------- */
6400 /* MAERMSG R*8 DSQRT I*4 MIN */
6403 /* DESCRIPTION/NOTES/LIMITATIONS : */
6404 /* ----------------------------------- */
6406 /* See VGAUSS to understand well the technique. */
6407 /* Actually SQRT (dpi^2) is integrated for i=1,nbdime */
6408 /* Calculation of the derivative is included in the code to avoid an additional */
6409 /* call of the routine. */
6411 /* The integrated function is strictly increasing, it */
6412 /* is not necessary to use a high degree for the GAUSS method GAUSS. */
6414 /* The degree of LEGENDRE polynom results from the degree of the */
6415 /* polynom to be integrated. It can vary from 4 to 40 (with step of 4). */
6417 /* The precision (relative) of integration is of order 1.D-8. */
6419 /* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
6421 /* Attention : the precision of the result is not controlled. */
6422 /* If you wish to control it, use MMCGLC1, taking into account that */
6423 /* the performance (in time) will be worse. */
6425 /* >=====================================================================
6428 /* ATTENTION : SAVE KGAR WGAUS and UROOT EVENTUALLY */
6430 /* INTEGER I1,I20 */
6431 /* PARAMETER (I1=1,I20=20) */
6433 /* Parameter adjustments */
6434 courbe_dim1 = *ndimax;
6435 courbe_offset = courbe_dim1 + 1;
6436 courbe -= courbe_offset;
6440 /* ****** General initialization ** */
6445 /* ****** Initialization of UROOT, WGAUS, NGAUS and KGAR ** */
6447 /* CALL MXVINIT(IERXV,'INTEGER',I1,KGAR,'INTEGER',I1,NGAUS */
6448 /* 1 ,'DOUBLE PRECISION',I20,UROOT,'DOUBLE PRECISION',I20,WGAUS) */
6449 /* IF (IERXV.GT.0) KGAR=0 */
6451 /* ****** Test the equity of limits ** */
6453 if (*tdebut == *tfinal) {
6458 /* ****** Test the dimension and the number of coefficients ** */
6460 if (*ndimen <= 0 || *ncoeff <= 0) {
6465 /* ****** Calculate the optimal degree ** */
6467 kk = *ncoeff / 4 + 1;
6470 /* ****** Return the coefficients for the integral (DEGRE=4*KK) */
6471 /* if KK <> KGAR. */
6474 mvgaus0_(&kk, uroot, wgaus, &ngaus, iercod);
6483 /* C1 => Point medium interval */
6484 /* C2 => 1/2 amplitude interval */
6486 c1 = (*tfinal + *tdebut) * .5;
6487 c2 = (*tfinal - *tdebut) * .5;
6489 /* ----------------------------------------------------------- */
6490 /* ****** Integration - Loop on GAUSS intervals ** */
6491 /* ----------------------------------------------------------- */
6496 for (jj = 1; jj <= i__1; ++jj) {
6498 /* ****** Integration taking the symmetry into account ** */
6500 tran = c2 * uroot[jj - 1];
6504 /* ****** Derivation on the dimension of the space ** */
6509 for (kk = 1; kk <= i__2; ++kk) {
6510 d1 = (*ncoeff - 1) * courbe[kk + *ncoeff * courbe_dim1];
6512 for (ii = *ncoeff - 1; ii >= 2; --ii) {
6513 dd = (ii - 1) * courbe[kk + ii * courbe_dim1];
6523 /* ****** Integration ** */
6525 som += wgaus[jj - 1] * c2 * (sqrt(der1) + sqrt(der2));
6527 /* ****** End of loop on GAUSS intervals ** */
6532 /* ****** Work ended ** */
6536 /* ****** It is forced IERCOD = 0 ** */
6540 /* ****** Final processing ** */
6544 /* ****** Save UROOT, WGAUS, NGAUS and KGAR ** */
6546 /* CALL MXVSAVE(IERXV,'INTEGER',I1,KGAR,'INTEGER',I1,NGAUS */
6547 /* 1 ,'DOUBLE PRECISION',I20,UROOT,'DOUBLE PRECISION',I20,WGAUS) */
6548 /* IF (IERXV.GT.0) KGAR=0 */
6550 /* ****** End of sub-program ** */
6553 AdvApp2Var_SysBase::maermsg_("MMLONCV", iercod, 7L);
6558 //=======================================================================
6559 //function : AdvApp2Var_MathBase::mmpobas_
6561 //=======================================================================
6562 int AdvApp2Var_MathBase::mmpobas_(doublereal *tparam,
6570 static integer c__2 = 2;
6571 static integer c__1 = 1;
6574 /* Initialized data */
6576 static doublereal moin11[2] = { -1.,1. };
6578 /* System generated locals */
6579 integer valbas_dim1, i__1;
6581 /* Local variables */
6582 static doublereal vjac[80], herm[24];
6583 static integer iord[2];
6584 static doublereal wval[4];
6585 static integer nwcof, iunit;
6586 static doublereal wpoly[7];
6587 static integer ii, jj, iorjac;
6588 static doublereal hermit[36] /* was [6][3][2] */;
6589 static integer kk1, kk2, kk3;
6590 static integer khe, ier;
6593 /* ***********************************************************************
6598 /* Position on the polynoms of base hermit-Jacobi */
6599 /* and their succesive derivatives */
6603 /* PUBLIC, POSITION, HERMIT, JACOBI */
6605 /* INPUT ARGUMENTS : */
6606 /* -------------------- */
6607 /* TPARAM : Parameter for which the position is found. */
6608 /* IORDRE : Orderof hermit-Jacobi (-1,0,1, ou 2) */
6609 /* NCOEFF : Number of coefficients of polynoms (Nb of value to calculate) */
6610 /* NDERIV : Number of derivative to calculate (0<= N <=3) */
6611 /* 0 -> Position simple on base functions */
6612 /* N -> Position on base functions and derivative */
6613 /* of order 1 to N */
6615 /* OUTPUT ARGUMENTS : */
6616 /* --------------------- */
6617 /* VALBAS (NCOEFF, 0:NDERIV) : calculated value */
6619 /* d vj(t) = VALBAS(J, I) */
6623 /* IERCOD : Error code */
6625 /* 1 : Incoherence of input arguments */
6627 /* COMMONS USED : */
6628 /* -------------- */
6631 /* REFERENCES CALLED : */
6632 /* ------------------- */
6635 /* DESCRIPTION/NOTES/LIMITATIONS : */
6636 /* ----------------------------------- */
6639 /* ***********************************************************************
6642 /* ***********************************************************************
6647 /* Parameter adjustments */
6648 valbas_dim1 = *ncoeff;
6653 /* ***********************************************************************
6655 /* INITIALIZATIONS */
6656 /* ***********************************************************************
6661 /* ***********************************************************************
6664 /* ***********************************************************************
6679 iorjac = (*iordre + 1) << 1;
6681 /* (1) Generic Calculations .... */
6683 /* (1.a) Calculation of hermit polynoms */
6686 mmherm1_(moin11, &c__2, iord, hermit, &ier);
6692 /* (1.b) Evaluation of hermit polynoms */
6695 iunit = *nderiv + 1;
6696 khe = (*iordre + 1) * iunit;
6701 for (ii = 0; ii <= i__1; ++ii) {
6702 mmdrvcb_(nderiv, &c__1, &iorjac, &hermit[(ii + 3) * 6 - 18],
6703 tparam, &herm[jj - 1], &ier);
6708 mmdrvcb_(nderiv, &c__1, &iorjac, &hermit[(ii + 6) * 6 - 18],
6709 tparam, &herm[jj + khe - 1], &ier);
6719 for (ii = 0; ii <= i__1; ++ii) {
6720 AdvApp2Var_MathBase::mmpocrb_(&c__1, &iorjac, &hermit[(ii + 3) * 6 - 18], &c__1,
6721 tparam, &herm[jj - 1]);
6723 AdvApp2Var_MathBase::mmpocrb_(&c__1, &iorjac, &hermit[(ii + 6) * 6 - 18], &c__1,
6724 tparam, &herm[jj + khe - 1]);
6729 /* (1.c) Evaluation of Jacobi polynoms */
6731 ii = *ncoeff - iorjac;
6733 mmpojac_(tparam, &iorjac, &ii, nderiv, vjac, &ier);
6738 /* (1.d) Evaluation of W(t) */
6742 nwcof = max(i__1,1);
6743 AdvApp2Var_SysBase::mvriraz_((integer *)&nwcof,
6750 } else if (*iordre == 1) {
6753 } else if (*iordre == 0) {
6757 mmdrvcb_(nderiv, &c__1, &nwcof, wpoly, tparam, wval, &ier);
6762 kk1 = *ncoeff - iorjac;
6766 /* (2) Evaluation of order 0 */
6770 for (ii = 1; ii <= i__1; ++ii) {
6771 valbas[ii] = herm[jj - 1];
6776 for (ii = 1; ii <= i__1; ++ii) {
6777 valbas[ii + iorjac] = wval[0] * vjac[ii - 1];
6780 /* (3) Evaluation of order 1 */
6785 for (ii = 1; ii <= i__1; ++ii) {
6786 valbas[ii + valbas_dim1] = herm[jj - 1];
6792 for (ii = 1; ii <= i__1; ++ii) {
6793 valbas[ii + iorjac + valbas_dim1] = wval[0] * vjac[ii + kk1 - 1]
6794 + wval[1] * vjac[ii - 1];
6798 /* (4) Evaluation of order 2 */
6803 for (ii = 1; ii <= i__1; ++ii) {
6804 valbas[ii + (valbas_dim1 << 1)] = herm[jj - 1];
6809 for (ii = 1; ii <= i__1; ++ii) {
6810 valbas[ii + iorjac + (valbas_dim1 << 1)] = wval[0] * vjac[ii +
6811 kk2 - 1] + wval[1] * 2 * vjac[ii + kk1 - 1] + wval[2] *
6816 /* (5) Evaluation of order 3 */
6821 for (ii = 1; ii <= i__1; ++ii) {
6822 valbas[ii + valbas_dim1 * 3] = herm[jj - 1];
6827 for (ii = 1; ii <= i__1; ++ii) {
6828 valbas[ii + iorjac + valbas_dim1 * 3] = wval[0] * vjac[ii + kk3 -
6829 1] + wval[1] * 3 * vjac[ii + kk2 - 1] + wval[2] * 3 *
6830 vjac[ii + kk1 - 1] + wval[3] * vjac[ii - 1];
6836 /* ***********************************************************************
6838 /* ERROR PROCESSING */
6839 /* ***********************************************************************
6849 /* ***********************************************************************
6851 /* RETURN CALLING PROGRAM */
6852 /* ***********************************************************************
6858 AdvApp2Var_SysBase::maermsg_("MMPOBAS", iercod, 7L);
6863 //=======================================================================
6864 //function : AdvApp2Var_MathBase::mmpocrb_
6866 //=======================================================================
6867 int AdvApp2Var_MathBase::mmpocrb_(integer *ndimax,
6875 /* System generated locals */
6876 integer courbe_dim1, courbe_offset, i__1, i__2;
6878 /* Local variables */
6879 static integer ncof2;
6880 static integer isize, nd, kcf, ncf;
6883 /* ***********************************************************************
6888 /* CALCULATE THE COORDINATES OF A POINT OF A CURVE OF GIVEN PARAMETER */
6889 /* TPARAM ( IN 2D, 3D OR MORE) */
6893 /* TOUS , MATH_ACCES :: COURBE&,PARAMETRE& , POSITIONNEMENT , &POINT
6896 /* INPUT ARGUMENTS : */
6897 /* ------------------ */
6898 /* NDIMAX : format / dimension of the curve */
6899 /* NCOEFF : Nb of coefficients of the curve */
6900 /* COURBE : Matrix of coefficients of the curve */
6901 /* NDIM : Dimension useful of the workspace */
6902 /* TPARAM : Value of the parameter where the point is calculated */
6904 /* OUTPUT ARGUMENTS : */
6905 /* ------------------- */
6906 /* PNTCRB : Coordinates of the calculated point */
6908 /* COMMONS USED : */
6909 /* ---------------- */
6913 /* REFERENCES CALLED : */
6914 /* ---------------------- */
6916 /* MIRAZ MVPSCR2 MVPSCR3 */
6918 /* DESCRIPTION/NOTES/LIMITATIONS : */
6919 /* ----------------------------------- */
6922 /* ***********************************************************************
6926 /* ***********************************************************************
6929 /* Parameter adjustments */
6930 courbe_dim1 = *ndimax;
6931 courbe_offset = courbe_dim1 + 1;
6932 courbe -= courbe_offset;
6937 AdvApp2Var_SysBase::miraz_((integer *)&isize,
6938 (char *)&pntcrb[1]);
6944 /* optimal processing 3d */
6946 if (*ndim == 3 && *ndimax == 3) {
6947 mvpscr3_(ncoeff, &courbe[courbe_offset], tparam, &pntcrb[1]);
6949 /* optimal processing 2d */
6951 } else if (*ndim == 2 && *ndimax == 2) {
6952 mvpscr2_(ncoeff, &courbe[courbe_offset], tparam, &pntcrb[1]);
6954 /* Any dimension - scheme of HORNER */
6956 } else if (*tparam == 0.) {
6958 for (nd = 1; nd <= i__1; ++nd) {
6959 pntcrb[nd] = courbe[nd + courbe_dim1];
6962 } else if (*tparam == 1.) {
6964 for (ncf = 1; ncf <= i__1; ++ncf) {
6966 for (nd = 1; nd <= i__2; ++nd) {
6967 pntcrb[nd] += courbe[nd + ncf * courbe_dim1];
6973 ncof2 = *ncoeff + 2;
6975 for (nd = 1; nd <= i__1; ++nd) {
6977 for (ncf = 2; ncf <= i__2; ++ncf) {
6979 pntcrb[nd] = (pntcrb[nd] + courbe[nd + kcf * courbe_dim1]) * *
6983 pntcrb[nd] += courbe[nd + courbe_dim1];
6992 //=======================================================================
6993 //function : AdvApp2Var_MathBase::mmmpocur_
6995 //=======================================================================
6996 int AdvApp2Var_MathBase::mmmpocur_(integer *ncofmx,
7004 /* System generated locals */
7005 integer courbe_dim1, courbe_offset, i__1;
7007 /* Local variables */
7008 static integer i__, nd;
7009 static doublereal fu;
7012 /* ***********************************************************************
7017 /* Position of a point on curve (ncofmx,ndim). */
7021 /* TOUS , AB_SPECIFI :: COURBE&,POLYNOME&,POSITIONNEMENT,&POINT */
7023 /* INPUT ARGUMENTS : */
7024 /* ------------------ */
7025 /* NCOFMX: Format / degree of the CURVE. */
7026 /* NDIM : Dimension of the space. */
7027 /* NDEG : Degree of the polynom. */
7028 /* COURBE: Coefficients of the curve. */
7029 /* TPARAM: Parameter on the curve */
7031 /* OUTPUT ARGUMENTS : */
7032 /* ------------------- */
7033 /* TABVAL(NDIM): The resulting point (or table of values) */
7035 /* COMMONS USED : */
7036 /* ---------------- */
7038 /* REFERENCES CALLED : */
7039 /* ----------------------- */
7041 /* DESCRIPTION/NOTES/LIMITATIONS : */
7042 /* ----------------------------------- */
7045 /* ***********************************************************************
7048 /* Parameter adjustments */
7050 courbe_dim1 = *ncofmx;
7051 courbe_offset = courbe_dim1 + 1;
7052 courbe -= courbe_offset;
7057 for (nd = 1; nd <= i__1; ++nd) {
7063 for (nd = 1; nd <= i__1; ++nd) {
7064 fu = courbe[*ndeg + nd * courbe_dim1];
7065 for (i__ = *ndeg - 1; i__ >= 1; --i__) {
7066 fu = fu * *tparam + courbe[i__ + nd * courbe_dim1];
7076 //=======================================================================
7077 //function : mmpojac_
7079 //=======================================================================
7080 int mmpojac_(doublereal *tparam,
7088 static integer c__2 = 2;
7090 /* Initialized data */
7092 static integer nbcof = -1;
7094 /* System generated locals */
7095 integer valjac_dim1, i__1, i__2;
7097 /* Local variables */
7098 static doublereal cofa, cofb, denom, tnorm[100];
7099 static integer ii, jj, kk1, kk2;
7100 static doublereal aux1, aux2;
7103 /* ***********************************************************************
7108 /* Positioning on Jacobi polynoms and their derivatives */
7109 /* successive by a recurrent algorithm */
7113 /* RESERVE, POSITIONING, JACOBI */
7115 /* INPUT ARGUMENTS : */
7116 /* -------------------- */
7117 /* TPARAM : Parameter for which positioning is done. */
7118 /* IORDRE : Order of hermit-?? (-1,0,1, or 2) */
7119 /* NCOEFF : Number of coeeficients of polynoms (Nb of value to */
7121 /* NDERIV : Number of derivative to calculate (0<= N <=3) */
7122 /* 0 -> Position simple on jacobi functions */
7123 /* N -> Position on jacobi functions and their */
7124 /* derivatives of order 1 to N. */
7126 /* OUTPUT ARGUMENTS : */
7127 /* --------------------- */
7128 /* VALJAC (NCOEFF, 0:NDERIV) : the calculated values */
7130 /* d vj(t) = VALJAC(J, I) */
7134 /* IERCOD : Error Code */
7136 /* 1 : Incoherence of input arguments */
7138 /* COMMONS USED : */
7139 /* ------------------ */
7142 /* REFERENCES CALLED : */
7143 /* --------------------- */
7146 /* DESCRIPTION/NOTES/LIMITATIONS : */
7147 /* ----------------------------------- */
7150 /* ***********************************************************************
7153 /* ***********************************************************************
7157 /* static varaibles */
7161 /* Parameter adjustments */
7162 valjac_dim1 = *ncoeff;
7167 /* ***********************************************************************
7169 /* INITIALISATIONS */
7170 /* ***********************************************************************
7175 /* ***********************************************************************
7178 /* ***********************************************************************
7184 if (*ncoeff > 100) {
7188 /* --- Calculation of norms */
7190 /* IF (NCOEFF.GT.NBCOF) THEN */
7192 for (ii = 1; ii <= i__1; ++ii) {
7196 for (jj = 1; jj <= i__2; ++jj) {
7197 aux2 = aux2 * (doublereal) (kk1 + *iordre + jj) / (doublereal) (
7200 i__2 = (*iordre << 1) + 1;
7201 tnorm[ii - 1] = sqrt(aux2 * (kk1 * 2. + (*iordre << 1) + 1) / pow__ii(&
7209 /* --- Trivial Positions ----- */
7212 aux1 = (doublereal) (*iordre + 1);
7213 valjac[2] = aux1 * *tparam;
7216 valjac[valjac_dim1 + 1] = 0.;
7217 valjac[valjac_dim1 + 2] = aux1;
7220 valjac[(valjac_dim1 << 1) + 1] = 0.;
7221 valjac[(valjac_dim1 << 1) + 2] = 0.;
7224 valjac[valjac_dim1 * 3 + 1] = 0.;
7225 valjac[valjac_dim1 * 3 + 2] = 0.;
7230 /* --- Positioning by recurrence */
7233 for (ii = 3; ii <= i__1; ++ii) {
7237 aux1 = (doublereal) (*iordre + kk2);
7239 cofa = aux2 * (aux2 + 1) * (aux2 + 2);
7240 cofb = (aux2 + 2) * -2. * aux1 * aux1;
7241 denom = kk1 * 2. * (kk2 + (*iordre << 1) + 1) * aux2;
7245 valjac[ii] = (cofa * *tparam * valjac[kk1] + cofb * valjac[kk2]) *
7249 valjac[ii + valjac_dim1] = (cofa * *tparam * valjac[kk1 +
7250 valjac_dim1] + cofa * valjac[kk1] + cofb * valjac[kk2 +
7251 valjac_dim1]) * denom;
7254 valjac[ii + (valjac_dim1 << 1)] = (cofa * *tparam * valjac[
7255 kk1 + (valjac_dim1 << 1)] + cofa * 2 * valjac[kk1 +
7256 valjac_dim1] + cofb * valjac[kk2 + (valjac_dim1 << 1)]
7261 valjac[ii + valjac_dim1 * 3] = (cofa * *tparam * valjac[kk1 +
7262 valjac_dim1 * 3] + cofa * 3 * valjac[kk1 + (
7263 valjac_dim1 << 1)] + cofb * valjac[kk2 + valjac_dim1 *
7269 /* ---> Normalization */
7272 for (ii = 1; ii <= i__1; ++ii) {
7274 for (jj = 0; jj <= i__2; ++jj) {
7275 valjac[ii + jj * valjac_dim1] = tnorm[ii - 1] * valjac[ii + jj *
7282 /* ***********************************************************************
7284 /* PROCESSING OF ERRORS */
7285 /* ***********************************************************************
7293 /* ***********************************************************************
7295 /* RETURN CALLING PROGRAM */
7296 /* ***********************************************************************
7302 AdvApp2Var_SysBase::maermsg_("MMPOJAC", iercod, 7L);
7307 //=======================================================================
7308 //function : AdvApp2Var_MathBase::mmposui_
7310 //=======================================================================
7311 int AdvApp2Var_MathBase::mmposui_(integer *dimmat,
7318 /* System generated locals */
7321 /* Local variables */
7322 static logical ldbg;
7323 static integer imin, jmin, i__, j, k;
7324 static logical trouve;
7326 /* ***********************************************************************
7331 /* FILL THE TABLE OF POSITIONING POSUIV WHICH ALLOWS TO */
7332 /* PARSE BY COLUMN THE INFERIOR TRIANGULAR PART OF THE */
7333 /* MATRIX IN FORM OF PROFILE */
7338 /* RESERVE, MATRIX, PROFILE */
7340 /* INPUT ARGUMENTS : */
7341 /* -------------------- */
7343 /* NISTOC: NUMBER OF COEFFICIENTS IN THE PROFILE */
7344 /* DIMMAT: NUMBER OF LINE OF THE SYMMETRIC SQUARE MATRIX */
7345 /* APOSIT: TABLE OF POSITIONING OF STORAGE TERMS */
7346 /* APOSIT(1,I) CONTAINS THE NUMBER OF TERMES-1 ON LINE
7347 /* I IN THE PROFILE OF THE MATRIX */
7348 /* APOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF DIAGONAL TERM
7352 /* OUTPUT ARGUMENTS : */
7353 /* --------------------- */
7354 /* POSUIV: POSUIV(K) (WHERE K IS THE INDEX OF STORAGE OF MAT(I,J)) */
7355 /* CONTAINS THE SMALLEST NUMBER IMIN>I OF THE LINE THAT */
7356 /* POSSESSES A TERM MAT(IMIN,J) THAT IS IN THE PROFILE. */
7357 /* IF THERE IS NO TERM MAT(IMIN,J) IN THE PROFILE THEN POSUIV(K)=-1 */
7360 /* COMMONS USED : */
7361 /* ------------------ */
7364 /* REFERENCES CALLED : */
7365 /* --------------------- */
7368 /* DESCRIPTION/NOTES/LIMITATIONS : */
7369 /* ----------------------------------- */
7372 /* ***********************************************************************
7375 /* ***********************************************************************
7380 /* ***********************************************************************
7382 /* INITIALIZATIONS */
7383 /* ***********************************************************************
7386 /* Parameter adjustments */
7391 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
7393 AdvApp2Var_SysBase::mgenmsg_("MMPOSUI", 7L);
7398 /* ***********************************************************************
7401 /* ***********************************************************************
7407 for (i__ = 1; i__ <= i__1; ++i__) {
7408 jmin = i__ - aposit[(i__ << 1) + 1];
7410 for (j = jmin; j <= i__2; ++j) {
7413 while(! trouve && imin <= *dimmat) {
7414 if (imin - aposit[(imin << 1) + 1] <= j) {
7420 k = aposit[(i__ << 1) + 2] - i__ + j;
7435 /* ***********************************************************************
7437 /* ERROR PROCESSING */
7438 /* ***********************************************************************
7444 /* ***********************************************************************
7446 /* RETURN CALLING PROGRAM */
7447 /* ***********************************************************************
7452 /* ___ DESALLOCATION, ... */
7454 AdvApp2Var_SysBase::maermsg_("MMPOSUI", iercod, 7L);
7456 AdvApp2Var_SysBase::mgsomsg_("MMPOSUI", 7L);
7461 //=======================================================================
7462 //function : AdvApp2Var_MathBase::mmresol_
7464 //=======================================================================
7465 int AdvApp2Var_MathBase::mmresol_(integer *hdimen,
7483 static integer c__100 = 100;
7485 /* System generated locals */
7488 /* Local variables */
7489 static logical ldbg;
7490 static doublereal mcho[100];
7491 static integer jmin, jmax, i__, j, k, l;
7492 static long int iofv1, iofv2, iofv3, iofv4;
7493 static doublereal v1[100], v2[100], v3[100], v4[100];
7494 static integer deblig, dimhch;
7495 static doublereal hchole[100];
7496 static long int iofmch, iofmam, iofhch;
7497 static doublereal matsym[100];
7503 /* ***********************************************************************
7508 /* SOLUTION OF THE SYSTEM */
7515 /* RESERVE, SOLUTION, SYSTEM, LAGRANGIAN */
7517 /* INPUT ARGUMENTS : */
7518 /* -------------------- */
7519 /* HDIMEN: NOMBER OF LINE (OR COLUMN) OF THE HESSIAN MATRIX */
7520 /* GDIMEN: NOMBER OF LINE OF THE MATRIX OF CONSTRAINTS */
7521 /* HNSTOC: NOMBErS OF TERMS IN THE PROFILE OF HESSIAN MATRIX
7523 /* GNSTOC: NOMBERS OF TERMS IN THE PROFILE OF THE MATRIX OF CONSTRAINTS */
7524 /* MNSTOC: NOMBERS OF TERMS IN THE PROFILE OF THE MATRIX M= G H t(G) */
7525 /* where H IS THE HESSIAN MATRIX AND G IS THE MATRIX OF CONSTRAINTS */
7526 /* MATSYH: TRIANGULAR INFERIOR PART OF THE HESSIAN MATRIX
7527 /* IN FORM OF PROFILE */
7528 /* MATSYG: MATRIX OF CONSTRAINTS IN FORM OF PROFILE */
7529 /* VECSYH: VECTOR OF THE SECOND MEMBER ASSOCIATED TO MATSYH */
7530 /* VECSYG: VECTOR OF THE SECOND MEMBER ASSOCIATED TO MATSYG */
7531 /* HPOSIT: TABLE OF POSITIONING OF THE HESSIAN MATRIX */
7532 /* HPOSIT(1,I) CONTAINS THE NUMBER OF TERMS -1 */
7533 /* WHICH ARE IN THE PROFILE AT LINE I */
7534 /* HPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF TERM */
7535 /* DIAGONAL OF THE MATRIX AT LINE I */
7536 /* HPOSUI: TABLE ALLOWING TO PARSE THE HESSIAN MATRIX BY COLUMN */
7537 /* IN FORM OF PROFILE */
7538 /* HPOSUI(K) CONTAINS THE NUMBER OF LINE IMIN FOLLOWING THE CURRENT LINE*/
7539 /* I WHERE H(I,J)=MATSYH(K) AS IT EXISTS IN THE */
7540 /* SAME COLUMN J A TERM IN THE PROFILE OF LINE IMIN */
7541 /* IF SUCH TERM DOES NOT EXIST IMIN=-1 */
7542 /* GPOSIT: TABLE OF POSITIONING OF THE MATRIX OF CONSTRAINTS */
7543 /* GPOSIT(1,I) CONTAINS THE NUMBER OF TERMS OF LINE I */
7544 /* WHICH ARE IN THE PROFILE */
7545 /* GPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF THE LAST TERM
7546 /* OF LINE I WHICH IS IN THE PROFILE */
7547 /* GPOSIT(3,I) CONTAINS THE NUMBER OF COLUMN CORRESPONDING */
7548 /* TO THE FIRST TERM OF LINE I WHICH IS IN THE PROFILE */
7549 /* MMPOSUI, MPOSIT: SAME STRUCTURE AS HPOSUI, BUT FOR MATRIX
7553 /* OUTPUT ARGUMENTS : */
7554 /* --------------------- */
7555 /* VECSOL: VECTOR SOLUTION V OF THE SYSTEM */
7556 /* IERCOD: ERROR CODE */
7558 /* COMMONS USED : */
7559 /* ------------------ */
7562 /* REFERENCES CALLED : */
7563 /* --------------------- */
7566 /* DESCRIPTION/NOTES/LIMITATIONS : */
7567 /* ----------------------------------- */
7569 /* ***********************************************************************
7572 /* ***********************************************************************
7575 /* ***********************************************************************
7577 /* INITIALISATIONS */
7578 /* ***********************************************************************
7581 /* Parameter adjustments */
7594 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
7596 AdvApp2Var_SysBase::mgenmsg_("MMRESOL", 7L);
7607 /* ***********************************************************************
7610 /* ***********************************************************************
7613 /* Dynamic allocation */
7615 AdvApp2Var_SysBase::macrar8_(hdimen, &c__100, v1, &iofv1, &ier);
7619 dimhch = hposit[(*hdimen << 1) + 2];
7620 AdvApp2Var_SysBase::macrar8_(&dimhch, &c__100, hchole, &iofhch, &ier);
7625 /* solution of system 1 H V1 = b */
7626 /* where H=MATSYH and b=VECSYH */
7628 mmchole_(hnstoc, hdimen, &matsyh[1], &hposit[3], &hposui[1], &hchole[
7633 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1], &vecsyh[
7634 1], &v1[iofv1], &ier);
7639 /* Case when there are constraints */
7643 /* Calculate the vector of the second member V2=G H(-1) b -c = G v1-c */
7644 /* of system of unknown Lagrangian vector MULTIP */
7645 /* where G=MATSYG */
7648 AdvApp2Var_SysBase::macrar8_(gdimen, &c__100, v2, &iofv2, &ier);
7652 AdvApp2Var_SysBase::macrar8_(hdimen, &c__100, v3, &iofv3, &ier);
7656 AdvApp2Var_SysBase::macrar8_(gdimen, &c__100, v4, &iofv4, &ier);
7660 AdvApp2Var_SysBase::macrar8_(mnstoc, &c__100, matsym, &iofmam, &ier);
7666 mmatvec_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v1[iofv1], &
7667 deblig, &v2[iofv2], &ier);
7672 for (i__ = 1; i__ <= i__1; ++i__) {
7673 v2[i__ + iofv2 - 1] -= vecsyg[i__];
7676 /* Calculate the matrix M= G H(-1) t(G) */
7677 /* RESOL DU SYST 2 : H qi = gi */
7678 /* where is a vector column of t(G) */
7680 /* then calculate G qi */
7681 /* then construct M in form of profile */
7686 for (i__ = 1; i__ <= i__1; ++i__) {
7687 AdvApp2Var_SysBase::mvriraz_((integer *)hdimen, (char *)&v1[iofv1]);
7688 AdvApp2Var_SysBase::mvriraz_((integer *)hdimen, (char *)&v3[iofv3]);
7689 AdvApp2Var_SysBase::mvriraz_((integer *)gdimen, (char *)&v4[iofv4]);
7690 jmin = gposit[i__ * 3 + 3];
7691 jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
7692 aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
7694 for (j = jmin; j <= i__2; ++j) {
7696 v1[j + iofv1 - 1] = matsyg[k];
7698 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1],
7699 &v1[iofv1], &v3[iofv3], &ier);
7705 mmatvec_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v3[
7706 iofv3], &deblig, &v4[iofv4], &ier);
7711 k = mposit[(i__ << 1) + 2];
7712 matsym[k + iofmam - 1] = v4[i__ + iofv4 - 1];
7713 while(mmposui[k] > 0) {
7715 k = mposit[(l << 1) + 2] - l + i__;
7716 matsym[k + iofmam - 1] = v4[l + iofv4 - 1];
7721 /* SOLVE SYST 3 M L = V2 */
7725 AdvApp2Var_SysBase::mvriraz_((integer *)gdimen, (char *)&v4[iofv4]);
7726 AdvApp2Var_SysBase::macrar8_(mnstoc, &c__100, mcho, &iofmch, &ier);
7730 mmchole_(mnstoc, gdimen, &matsym[iofmam], &mposit[3], &mmposui[1], &
7731 mcho[iofmch], &ier);
7735 mmrslss_(mnstoc, gdimen, &mcho[iofmch], &mposit[3], &mmposui[1], &v2[
7736 iofv2], &v4[iofv4], &ier);
7742 /* CALCULATE THE VECTOR OF THE SECOND MEMBER OF THE SYSTEM Hx = b - t(G) L
7746 AdvApp2Var_SysBase::mvriraz_((integer *)hdimen, (char *)&v1[iofv1]);
7747 mmtmave_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v4[iofv4], &
7753 for (i__ = 1; i__ <= i__1; ++i__) {
7754 v1[i__ + iofv1 - 1] = vecsyh[i__] - v1[i__ + iofv1 - 1];
7757 /* RESOL SYST 4 Hx = b - t(G) L */
7760 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1], &v1[
7761 iofv1], &vecsol[1], &ier);
7767 for (i__ = 1; i__ <= i__1; ++i__) {
7768 vecsol[i__] = v1[i__ + iofv1 - 1];
7774 /* ***********************************************************************
7776 /* PROCESSING OF ERRORS */
7777 /* ***********************************************************************
7786 AdvApp2Var_SysBase::mswrdbg_("MMRESOL : PROBLEM WITH DIMMAT", 30L);
7789 /* ***********************************************************************
7791 /* RETURN CALLING PROGRAM */
7792 /* ***********************************************************************
7797 /* ___ DESALLOCATION, ... */
7798 AdvApp2Var_SysBase::macrdr8_(hdimen, &c__100, v1, &iofv1, &ier);
7799 if (*iercod == 0 && ier > 0) {
7802 AdvApp2Var_SysBase::macrdr8_(&dimhch, &c__100, hchole, &iofhch, &ier);
7803 if (*iercod == 0 && ier > 0) {
7806 AdvApp2Var_SysBase::macrdr8_(gdimen, &c__100, v2, &iofv2, &ier);
7807 if (*iercod == 0 && ier > 0) {
7810 AdvApp2Var_SysBase::macrdr8_(hdimen, &c__100, v3, &iofv3, &ier);
7811 if (*iercod == 0 && ier > 0) {
7814 AdvApp2Var_SysBase::macrdr8_(gdimen, &c__100, v4, &iofv4, &ier);
7815 if (*iercod == 0 && ier > 0) {
7818 AdvApp2Var_SysBase::macrdr8_(mnstoc, &c__100, matsym, &iofmam, &ier);
7819 if (*iercod == 0 && ier > 0) {
7822 AdvApp2Var_SysBase::macrdr8_(mnstoc, &c__100, mcho, &iofmch, &ier);
7823 if (*iercod == 0 && ier > 0) {
7827 AdvApp2Var_SysBase::maermsg_("MMRESOL", iercod, 7L);
7829 AdvApp2Var_SysBase::mgsomsg_("MMRESOL", 7L);
7834 //=======================================================================
7835 //function : mmrslss_
7837 //=======================================================================
7838 int mmrslss_(integer *,//mxcoef,
7843 doublereal *mscnmbr,
7847 /* System generated locals */
7850 /* Local variables */
7851 static logical ldbg;
7852 static integer i__, j;
7853 static doublereal somme;
7854 static integer pointe, ptcour;
7856 /* ***********************************************************************
7861 /* Solves linear system SS x = b where S is a */
7862 /* triangular lower matrix given in form of profile */
7866 /* RESERVE, MATRICE_PROFILE, RESOLUTION, CHOLESKI */
7868 /* INPUT ARGUMENTS : */
7869 /* -------------------- */
7870 /* MXCOEF : Maximum number of non-null coefficient in the matrix */
7871 /* DIMENS : Dimension of the matrix */
7872 /* SMATRI(MXCOEF) : Values of coefficients of the matrix */
7873 /* SPOSIT(2,DIMENS): */
7874 /* SPOSIT(1,*) : Distance diagonal-extremity of the line */
7875 /* SPOSIT(2,*) : Position of diagonal terms in AMATRI */
7876 /* POSUIV(MXCOEF): first line inferior not out of profile */
7877 /* MSCNMBR(DIMENS): Vector second member of the equation */
7879 /* OUTPUT ARGUMENTS : */
7880 /* --------------------- */
7881 /* SOLUTI(NDIMEN) : Result vector */
7882 /* IERCOD : Error code 0 : ok */
7884 /* COMMONS USED : */
7885 /* ------------------ */
7888 /* REFERENCES CALLED : */
7889 /* --------------------- */
7892 /* DESCRIPTION/NOTES/LIMITATIONS : */
7893 /* ----------------------------------- */
7895 /* SS is the decomposition of choleski of a symmetric matrix */
7896 /* defined postive, that can result from routine MMCHOLE. */
7898 /* For a full matrix it is possible to use MRSLMSC */
7900 /* LEVEL OF DEBUG = 4 */
7902 /* ***********************************************************************
7905 /* ***********************************************************************
7910 /* ***********************************************************************
7912 /* INITIALISATIONS */
7913 /* ***********************************************************************
7916 /* Parameter adjustments */
7924 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
7926 AdvApp2Var_SysBase::mgenmsg_("MMRSLSS", 7L);
7930 /* ***********************************************************************
7933 /* ***********************************************************************
7936 /* ----- Solution of Sw = b */
7939 for (i__ = 1; i__ <= i__1; ++i__) {
7941 pointe = sposit[(i__ << 1) + 2];
7944 for (j = i__ - sposit[(i__ << 1) + 1]; j <= i__2; ++j) {
7945 somme += smatri[pointe - (i__ - j)] * soluti[j];
7948 soluti[i__] = (mscnmbr[i__] - somme) / smatri[pointe];
7951 /* ----- Solution of S u = w */
7953 for (i__ = *dimens; i__ >= 1; --i__) {
7955 pointe = sposit[(i__ << 1) + 2];
7959 ptcour = sposit[(j << 1) + 2] - (j - i__);
7960 somme += smatri[ptcour] * soluti[j];
7964 soluti[i__] = (soluti[i__] - somme) / smatri[pointe];
7969 /* ***********************************************************************
7971 /* ERROR PROCESSING */
7972 /* ***********************************************************************
7976 /* ***********************************************************************
7978 /* RETURN PROGRAM CALLING */
7979 /* ***********************************************************************
7984 AdvApp2Var_SysBase::maermsg_("MMRSLSS", iercod, 7L);
7986 AdvApp2Var_SysBase::mgsomsg_("MMRSLSS", 7L);
7991 //=======================================================================
7992 //function : mmrslw_
7994 //=======================================================================
7995 int mmrslw_(integer *normax,
8003 /* System generated locals */
8004 integer abmatr_dim1, abmatr_offset, xmatri_dim1, xmatri_offset, i__1,
8008 /* Local variables */
8009 static integer kpiv;
8010 static doublereal pivot;
8011 static integer ii, jj, kk;
8012 static doublereal akj;
8015 /* **********************************************************************
8020 /* Solution of a linear system A.x = B of N equations to N */
8021 /* unknown by Gauss method (partial pivot) or : */
8022 /* A is matrix NORDRE * NORDRE, */
8023 /* B is matrix NORDRE (lines) * NDIMEN (columns), */
8024 /* x is matrix NORDRE (lines) * NDIMEN (columns). */
8025 /* In this program, A and B are stored in matrix ABMATR */
8026 /* the lines and columns which of were inverted. ABMATR(k,j) is */
8027 /* term A(j,k) if k <= NORDRE, B(j,k-NORDRE) otherwise (see example). */
8031 /* TOUS, MATH_ACCES::EQUATION&, MATRICE&, RESOLUTION, GAUSS, &SOLUTION */
8033 /* INPUT ARGUMENTS : */
8034 /* ------------------ */
8035 /* NORMAX : Max size of the first index of XMATRI. This argument */
8036 /* serves only for the declaration of dimension of XMATRI and should be */
8037 /* above or equal to NORDRE. */
8038 /* NORDRE : Order of the matrix i.e. number of equations and */
8039 /* unknown quantities of the linear system to be solved. */
8040 /* NDIMEN : Number of the second member. */
8041 /* EPSPIV : Minimal value of a pivot. If during the calculation */
8042 /* the absolute value of the pivot is below EPSPIV, the */
8043 /* system of equations is declared singular. EPSPIV should */
8044 /* be a "small" real. */
8046 /* ABMATR(NORDRE+NDIMEN,NORDRE) : Auxiliary matrix containing */
8047 /* matrix A and matrix B. */
8049 /* OUTPUT ARGUMENTS : */
8050 /* ------------------- */
8051 /* XMATRI : Matrix containing NORDRE*NDIMEN solutions. */
8052 /* IERCOD=0 shows that all solutions are calculated. */
8053 /* IERCOD=1 shows that the matrix is of lower rank than NORDRE */
8054 /* (the system is singular). */
8056 /* COMMONS USED : */
8057 /* ---------------- */
8059 /* REFERENCES CALLED : */
8060 /* ----------------------- */
8062 /* DESCRIPTION/NOTES/LIMITATIONS : */
8063 /* ----------------------------------- */
8064 /* ATTENTION : the indices of line and column are inverted */
8065 /* compared to usual indices. */
8067 /* a1*x + b1*y = c1 */
8068 /* a2*x + b2*y = c2 */
8069 /* should be represented by matrix ABMATR : */
8071 /* ABMATR(1,1) = a1 ABMATR(1,2) = a2 */
8072 /* ABMATR(2,1) = b1 ABMATR(2,2) = b2 */
8073 /* ABMATR(3,1) = c1 ABMATR(3,2) = c2 */
8075 /* To solve this system, it is necessary to set : */
8077 /* NORDRE = 2 (there are 2 equations with 2 unknown values), */
8078 /* NDIMEN = 1 (there is only one second member), */
8079 /* any NORMAX can be taken >= NORDRE. */
8081 /* To use this routine, it is recommended to use one of */
8082 /* interfaces : MMRSLWI or MMMRSLWD. */
8084 /* **********************************************************************
8087 /* Name of the routine */
8089 /* INTEGER IBB,MNFNDEB */
8092 /* IF (IBB.GE.2) CALL MGENMSG(NOMPR) */
8093 /* Parameter adjustments */
8094 xmatri_dim1 = *normax;
8095 xmatri_offset = xmatri_dim1 + 1;
8096 xmatri -= xmatri_offset;
8097 abmatr_dim1 = *nordre + *ndimen;
8098 abmatr_offset = abmatr_dim1 + 1;
8099 abmatr -= abmatr_offset;
8104 /* *********************************************************************
8106 /* Triangulation of matrix ABMATR. */
8107 /* *********************************************************************
8111 for (kk = 1; kk <= i__1; ++kk) {
8113 /* ---------- Find max pivot in column KK. ------------
8119 for (jj = kk; jj <= i__2; ++jj) {
8120 akj = (d__1 = abmatr[kk + jj * abmatr_dim1], abs(d__1));
8131 /* --------- Swapping of line KPIV with line KK. ------
8135 i__2 = *nordre + *ndimen;
8136 for (jj = kk; jj <= i__2; ++jj) {
8137 akj = abmatr[jj + kk * abmatr_dim1];
8138 abmatr[jj + kk * abmatr_dim1] = abmatr[jj + kpiv *
8140 abmatr[jj + kpiv * abmatr_dim1] = akj;
8145 /* ---------- Removal and triangularization. -----------
8148 pivot = -abmatr[kk + kk * abmatr_dim1];
8150 for (ii = kk + 1; ii <= i__2; ++ii) {
8151 akj = abmatr[kk + ii * abmatr_dim1] / pivot;
8152 i__3 = *nordre + *ndimen;
8153 for (jj = kk + 1; jj <= i__3; ++jj) {
8154 abmatr[jj + ii * abmatr_dim1] += akj * abmatr[jj + kk *
8165 /* *********************************************************************
8167 /* Solution of the system of triangular equations. */
8168 /* Matrix ABMATR(NORDRE+JJ,II), contains second members */
8169 /* of the system for 1<=j<=NDIMEN and 1<=i<=NORDRE. */
8170 /* *********************************************************************
8174 /* ---------------- Calculation of solutions by ascending. -----------------
8177 for (kk = *nordre; kk >= 1; --kk) {
8178 pivot = abmatr[kk + kk * abmatr_dim1];
8180 for (ii = 1; ii <= i__1; ++ii) {
8181 akj = abmatr[ii + *nordre + kk * abmatr_dim1];
8183 for (jj = kk + 1; jj <= i__2; ++jj) {
8184 akj -= abmatr[jj + kk * abmatr_dim1] * xmatri[jj + ii *
8188 xmatri[kk + ii * xmatri_dim1] = akj / pivot;
8195 /* ------If the absolute value of a pivot is smaller than --------
8196 /* ---------- EPSPIV: return the code of error. ------------
8206 AdvApp2Var_SysBase::maermsg_("MMRSLW ", iercod, 7L);
8208 /* IF (IBB.GE.2) CALL MGSOMSG(NOMPR) */
8212 //=======================================================================
8213 //function : AdvApp2Var_MathBase::mmmrslwd_
8215 //=======================================================================
8216 int AdvApp2Var_MathBase::mmmrslwd_(integer *normax,
8227 /* System generated locals */
8228 integer amat_dim1, amat_offset, bmat_dim1, bmat_offset, xmat_dim1,
8229 xmat_offset, aaux_dim1, aaux_offset, i__1, i__2;
8231 /* Local variables */
8232 static integer i__, j;
8235 /* IMPLICIT DOUBLE PRECISION (A-H,O-Z) */
8236 /* IMPLICIT INTEGER (I-N) */
8239 /* **********************************************************************
8244 /* Solution of a linear system by Gauss method where */
8245 /* the second member is a table of vectors. Method of partial pivot. */
8249 /* ALL, MATH_ACCES :: */
8250 /* SYSTEME&,EQUATION&, RESOLUTION,GAUSS ,&VECTEUR */
8252 /* INPUT ARGUMENTS : */
8253 /* ------------------ */
8254 /* NORMAX : Max. Dimension of AMAT. */
8255 /* NORDRE : Order of the matrix. */
8256 /* NDIM : Number of columns of BMAT and XMAT. */
8257 /* AMAT(NORMAX,NORDRE) : The processed matrix. */
8258 /* BMAT(NORMAX,NDIM) : The matrix of second member. */
8259 /* XMAT(NORMAX,NDIM) : The matrix of solutions. */
8260 /* EPSPIV : Min value of a pivot. */
8262 /* OUTPUT ARGUMENTS : */
8263 /* ------------------- */
8264 /* AAUX(NORDRE+NDIM,NORDRE) : Auxiliary matrix. */
8265 /* XMAT(NORMAX,NDIM) : Matrix of solutions. */
8266 /* IERCOD=0 shows that solutions in XMAT are valid. */
8267 /* IERCOD=1 shows that matrix AMAT is of lower rank than NORDRE. */
8269 /* COMMONS USED : */
8270 /* ---------------- */
8274 /* REFERENCES CALLED : */
8275 /* ---------------------- */
8277 /* MAERMSG MGENMSG MGSOMSG */
8278 /* MMRSLW I*4 MNFNDEB */
8280 /* DESCRIPTION/NOTES/LIMITATIONS : */
8281 /* ----------------------------------- */
8282 /* ATTENTION : lines and columns are located in usual order : */
8283 /* 1st index = index line */
8284 /* 2nd index = index column */
8285 /* Example, the system : */
8286 /* a1*x + b1*y = c1 */
8287 /* a2*x + b2*y = c2 */
8288 /* is represented by matrix AMAT : */
8290 /* AMAT(1,1) = a1 AMAT(2,1) = a2 */
8291 /* AMAT(1,2) = b1 AMAT(2,2) = b2 */
8293 /* The first index is the index of line, the second index */
8294 /* is the index of columns (Compare with MMRSLWI which is faster). */
8297 /* **********************************************************************
8300 /* Name of the routine */
8302 /* Parameter adjustments */
8303 amat_dim1 = *normax;
8304 amat_offset = amat_dim1 + 1;
8305 amat -= amat_offset;
8306 xmat_dim1 = *normax;
8307 xmat_offset = xmat_dim1 + 1;
8308 xmat -= xmat_offset;
8309 aaux_dim1 = *nordre + *ndim;
8310 aaux_offset = aaux_dim1 + 1;
8311 aaux -= aaux_offset;
8312 bmat_dim1 = *normax;
8313 bmat_offset = bmat_dim1 + 1;
8314 bmat -= bmat_offset;
8317 ibb = AdvApp2Var_SysBase::mnfndeb_();
8319 AdvApp2Var_SysBase::mgenmsg_("MMMRSLW", 7L);
8322 /* Initialization of the auxiliary matrix. */
8325 for (i__ = 1; i__ <= i__1; ++i__) {
8327 for (j = 1; j <= i__2; ++j) {
8328 aaux[j + i__ * aaux_dim1] = amat[i__ + j * amat_dim1];
8334 /* Second member. */
8337 for (i__ = 1; i__ <= i__1; ++i__) {
8339 for (j = 1; j <= i__2; ++j) {
8340 aaux[j + *nordre + i__ * aaux_dim1] = bmat[i__ + j * bmat_dim1];
8346 /* Solution of the system of equations. */
8348 mmrslw_(normax, nordre, ndim, epspiv, &aaux[aaux_offset], &xmat[
8349 xmat_offset], iercod);
8353 AdvApp2Var_SysBase::maermsg_("MMMRSLW", iercod, 7L);
8356 AdvApp2Var_SysBase::mgsomsg_("MMMRSLW", 7L);
8361 //=======================================================================
8362 //function : AdvApp2Var_MathBase::mmrtptt_
8364 //=======================================================================
8365 int AdvApp2Var_MathBase::mmrtptt_(integer *ndglgd,
8369 static integer ideb, nmod2, nsur2, ilong, ibb;
8372 /* **********************************************************************
8377 /* Extracts from Common LDGRTL the STRICTLY positive roots of the */
8378 /* Legendre polynom of degree NDGLGD, for 2 <= NDGLGD <= 61. */
8382 /* TOUS, AB_SPECIFI::COMMON&, EXTRACTION, &RACINE, &LEGENDRE. */
8384 /* INPUT ARGUMENTS : */
8385 /* ------------------ */
8386 /* NDGLGD : Mathematic degree of Legendre polynom. */
8387 /* This degree should be above or equal to 2 and */
8388 /* below or equal to 61. */
8390 /* OUTPUT ARGUMENTS : */
8391 /* ------------------- */
8392 /* RTLEGD : The table of strictly positive roots of */
8393 /* Legendre polynom of degree NDGLGD. */
8395 /* COMMONS USED : */
8396 /* ---------------- */
8398 /* REFERENCES CALLED : */
8399 /* ----------------------- */
8401 /* DESCRIPTION/NOTES/LIMITATIONS : */
8402 /* ----------------------------------- */
8403 /* ATTENTION: the condition on NDEGRE ( 2 <= NDEGRE <= 61) is not */
8404 /* tested. The caller should make the test. */
8407 /* **********************************************************************
8409 /* Nome of the routine */
8412 /* Common MLGDRTL: */
8413 /* This common includes POSITIVE roots of Legendre polynoms */
8414 /* AND the weight of Gauss quadrature formulas on all */
8415 /* POSITIVE roots of Legendre polynoms. */
8418 /* ***********************************************************************
8423 /* The common of Legendre roots. */
8429 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
8430 /* ----------------------------------- */
8433 /* ***********************************************************************
8439 /* ROOTAB : Table of all rotts of Legendre polynoms */
8440 /* between [0,1]. They are ranked for degrees increasing from 2 to 61. */
8441 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
8442 /* The address is the same. */
8443 /* HI0TAB : Table of Legendre interpolators for root x=0 */
8444 /* the polynoms of UNEVEN degree. */
8445 /* RTLTB0 : Table of Li(uk) where uk are roots of a */
8446 /* Legendre polynom of EVEN degree. */
8447 /* RTLTB1 : Table of Li(uk) where uk are roots of a */
8448 /* Legendre polynom of UNEVEN degree. */
8451 /************************************************************************
8453 /* Parameter adjustments */
8457 ibb = AdvApp2Var_SysBase::mnfndeb_();
8459 AdvApp2Var_SysBase::mgenmsg_("MMRTPTT", 7L);
8465 nsur2 = *ndglgd / 2;
8466 nmod2 = *ndglgd % 2;
8469 ideb = nsur2 * (nsur2 - 1) / 2 + 1;
8470 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
8471 (char *)&mlgdrtl_.rootab[ideb + nmod2 * 465 - 1],
8472 (char *)&rtlegd[1]);
8474 /* ----------------------------- The end --------------------------------
8479 AdvApp2Var_SysBase::mgsomsg_("MMRTPTT", 7L);
8484 //=======================================================================
8485 //function : AdvApp2Var_MathBase::mmsrre2_
8487 //=======================================================================
8488 int AdvApp2Var_MathBase::mmsrre2_(doublereal *tparam,
8496 /* System generated locals */
8499 /* Local variables */
8500 static integer ideb, ifin, imil, ibb;
8502 /* ***********************************************************************
8508 /* Find the interval corresponding to a valueb given in */
8509 /* increasing order of real numbers with double precision. */
8513 /* TOUS,MATH_ACCES::TABLEAU&,POINT&,CORRESPONDANCE,&RANG */
8515 /* INPUT ARGUMENTS : */
8516 /* ------------------ */
8518 /* TPARAM : Value to be tested. */
8519 /* NBRVAL : Size of TABLEV */
8520 /* TABLEV : Table of reals. */
8521 /* EPSIL : Epsilon of precision */
8523 /* OUTPUT ARGUMENTS : */
8524 /* ------------------- */
8526 /* NUMINT : Number of the interval (between 1 and NBRVAL-1). */
8527 /* ITYPEN : = 0 TPARAM is inside the interval NUMINT */
8528 /* = 1 : TPARAM corresponds to the lower limit of */
8529 /* the provided interval. */
8530 /* = 2 : TPARAM corresponds to the upper limit of */
8531 /* the provided interval. */
8533 /* IERCOD : Error code. */
8535 /* = 1 : TABLEV does not contain enough elements. */
8536 /* = 2 : TPARAM out of limits of TABLEV. */
8538 /* COMMONS USED : */
8539 /* ---------------- */
8541 /* REFERENCES CALLED : */
8542 /* ------------------- */
8544 /* DESCRIPTION/NOTES/LIMITATIONS : */
8545 /* --------------------------------- */
8546 /* There are NBRVAL values in TABLEV which stands for NBRVAL-1 intervals. */
8547 /* One searches the interval containing TPARAM by */
8548 /* dichotomy. Complexity of the algorithm : Log(n)/Log(2).(RBD). */
8550 /* ***********************************************************************
8554 /* Initialisations */
8556 /* Parameter adjustments */
8560 ibb = AdvApp2Var_SysBase::mnfndeb_();
8562 AdvApp2Var_SysBase::mgenmsg_("MMSRRE2", 7L);
8571 /* TABLEV should contain at least two values */
8578 /* TPARAM should be between extreme limits of TABLEV. */
8580 if (*tparam < tablev[1] || *tparam > tablev[*nbrval]) {
8585 /* ----------------------- SEARCH OF THE INTERVAL --------------------
8590 /* Test end of loop (found). */
8592 if (ideb + 1 == ifin) {
8597 /* Find by dichotomy on increasing values of TABLEV. */
8599 imil = (ideb + ifin) / 2;
8600 if (*tparam >= tablev[ideb] && *tparam <= tablev[imil]) {
8608 /* -------------- TEST IF TPARAM IS NOT A VALUE ---------
8609 /* ------------------------OF TABLEV UP TO EPSIL ----------------------
8613 if ((d__1 = *tparam - tablev[ideb], abs(d__1)) < *epsil) {
8617 if ((d__1 = *tparam - tablev[ifin], abs(d__1)) < *epsil) {
8622 /* --------------------------- THE END ----------------------------------
8627 AdvApp2Var_SysBase::maermsg_("MMSRRE2", iercod, 7L);
8630 AdvApp2Var_SysBase::mgsomsg_("MMSRRE2", 7L);
8635 //=======================================================================
8636 //function : mmtmave_
8638 //=======================================================================
8639 int mmtmave_(integer *nligne,
8649 /* System generated locals */
8652 /* Local variables */
8653 static logical ldbg;
8654 static integer imin, imax, i__, j, k;
8655 static doublereal somme;
8659 /* ***********************************************************************
8665 /* CREATES PRODUCT G V */
8666 /* WHERE THE MATRIX IS IN FORM OF PROFILE */
8670 /* RESERVE, PRODUCT, MATRIX, PROFILE, VECTOR */
8672 /* INPUT ARGUMENTS : */
8673 /* -------------------- */
8674 /* NLIGNE : NUMBER OF LINE OF THE MATRIX */
8675 /* NCOLON : NOMBER OF COLUMN OF THE MATRIX */
8676 /* GPOSIT: TABLE OF POSITIONING OF TERMS OF STORAGE */
8677 /* GPOSIT(1,I) CONTAINS THE NUMBER of TERMS-1 ON LINE
8678 I IN THE PROFILE OF THE MATRIX */
8679 /* GPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF THE DIAGONAL TERM
8681 /* GPOSIT(3,I) CONTAINS THE INDEX COLUMN OF THE FIRST TERM OF
8682 /* PROFILE OF LINE I */
8683 /* GNSTOC : NOMBER OF TERM IN THE PROFILE OF GMATRI */
8684 /* GMATRI : MATRIX OF CONSTRAINTS IN FORM OF PROFILE */
8685 /* VECIN : INPUT VECTOR */
8687 /* OUTPUT ARGUMENTS : */
8688 /* --------------------- */
8689 /* VECOUT : VECTOR PRODUCT */
8690 /* IERCOD : ERROR CODE */
8693 /* COMMONS USED : */
8694 /* ------------------ */
8697 /* REFERENCES CALLED : */
8698 /* --------------------- */
8701 /* DESCRIPTION/NOTES/LIMITATIONS : */
8702 /* ----------------------------------- */
8704 /* ***********************************************************************
8707 /* ***********************************************************************
8712 /* ***********************************************************************
8714 /* INITIALISATIONS */
8715 /* ***********************************************************************
8718 /* Parameter adjustments */
8725 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
8727 AdvApp2Var_SysBase::mgenmsg_("MMTMAVE", 7L);
8731 /* ***********************************************************************
8734 /* ***********************************************************************
8740 for (i__ = 1; i__ <= i__1; ++i__) {
8743 for (j = 1; j <= i__2; ++j) {
8744 imin = gposit[j * 3 + 3];
8745 imax = gposit[j * 3 + 1] + gposit[j * 3 + 3] - 1;
8746 aux = gposit[j * 3 + 2] - gposit[j * 3 + 1] - imin + 1;
8747 if (imin <= i__ && i__ <= imax) {
8749 somme += gmatri[k] * vecin[j];
8752 vecout[i__] = somme;
8761 /* ***********************************************************************
8763 /* ERROR PROCESSING */
8764 /* ***********************************************************************
8768 /* ***********************************************************************
8770 /* RETURN CALLING PROGRAM */
8771 /* ***********************************************************************
8776 /* ___ DESALLOCATION, ... */
8778 AdvApp2Var_SysBase::maermsg_("MMTMAVE", iercod, 7L);
8780 AdvApp2Var_SysBase::mgsomsg_("MMTMAVE", 7L);
8785 //=======================================================================
8786 //function : mmtrpj0_
8788 //=======================================================================
8789 int mmtrpj0_(integer *ncofmx,
8799 /* System generated locals */
8800 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
8803 /* Local variables */
8804 static integer ncut, i__;
8805 static doublereal bidon, error;
8809 /* ***********************************************************************
8814 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
8815 /* Legendre with a given precision. */
8819 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
8821 /* INPUT ARGUMENTS : */
8822 /* ------------------ */
8823 /* NCOFMX : Max Nb of coeff. of the curve (dimensioning). */
8824 /* NDIMEN : Dimension of the space. */
8825 /* NCOEFF : Degree +1 of the polynom. */
8826 /* EPSI3D : Precision required for the approximation. */
8827 /* CRVLGD : The curve the degree which of it is required to lower. */
8829 /* OUTPUT ARGUMENTS : */
8830 /* ------------------- */
8831 /* EPSTRC : Precision of the approximation. */
8832 /* NCFNEW : Degree +1 of the resulting polynom. */
8834 /* COMMONS USED : */
8835 /* ---------------- */
8837 /* REFERENCES CALLED : */
8838 /* ----------------------- */
8840 /* DESCRIPTION/NOTES/LIMITATIONS : */
8841 /* ----------------------------------- */
8843 /* ***********************************************************************
8847 /* ------- Minimum degree that can be attained : Stop at 1 (RBD) ---------
8850 /* Parameter adjustments */
8852 crvlgd_dim1 = *ncofmx;
8853 crvlgd_offset = crvlgd_dim1 + 1;
8854 crvlgd -= crvlgd_offset;
8858 /* ------------------- Init for error calculation -----------------------
8861 for (i__ = 1; i__ <= i__1; ++i__) {
8868 /* Cutting of coefficients. */
8871 /* ------ Loop on the series of Legendre :NCOEFF --> 2 (RBD) -----------
8874 for (i__ = *ncoeff; i__ >= i__1; --i__) {
8875 /* Factor of renormalization. */
8876 bidon = ((i__ - 1) * 2. + 1.) / 2.;
8877 bidon = sqrt(bidon);
8879 for (nd = 1; nd <= i__2; ++nd) {
8880 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], abs(d__1)) *
8884 /* Cutting is stopped if the norm becomes too great. */
8885 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
8886 if (error > *epsi3d) {
8891 /* --- Max error cumulee when the I-th coeff is removed. */
8898 /* --------------------------------- End --------------------------------
8905 //=======================================================================
8906 //function : mmtrpj2_
8908 //=======================================================================
8909 int mmtrpj2_(integer *ncofmx,
8919 /* Initialized data */
8921 static doublereal xmaxj[57] = { .9682458365518542212948163499456,
8922 .986013297183269340427888048593603,
8923 1.07810420343739860362585159028115,
8924 1.17325804490920057010925920756025,
8925 1.26476561266905634732910520370741,
8926 1.35169950227289626684434056681946,
8927 1.43424378958284137759129885012494,
8928 1.51281316274895465689402798226634,
8929 1.5878364329591908800533936587012,
8930 1.65970112228228167018443636171226,
8931 1.72874345388622461848433443013543,
8932 1.7952515611463877544077632304216,
8933 1.85947199025328260370244491818047,
8934 1.92161634324190018916351663207101,
8935 1.98186713586472025397859895825157,
8936 2.04038269834980146276967984252188,
8937 2.09730119173852573441223706382076,
8938 2.15274387655763462685970799663412,
8939 2.20681777186342079455059961912859,
8940 2.25961782459354604684402726624239,
8941 2.31122868752403808176824020121524,
8942 2.36172618435386566570998793688131,
8943 2.41117852396114589446497298177554,
8944 2.45964731268663657873849811095449,
8945 2.50718840313973523778244737914028,
8946 2.55385260994795361951813645784034,
8947 2.59968631659221867834697883938297,
8948 2.64473199258285846332860663371298,
8949 2.68902863641518586789566216064557,
8950 2.73261215675199397407027673053895,
8951 2.77551570192374483822124304745691,
8952 2.8177699459714315371037628127545,
8953 2.85940333797200948896046563785957,
8954 2.90044232019793636101516293333324,
8955 2.94091151970640874812265419871976,
8956 2.98083391718088702956696303389061,
8957 3.02023099621926980436221568258656,
8958 3.05912287574998661724731962377847,
8959 3.09752842783622025614245706196447,
8960 3.13546538278134559341444834866301,
8961 3.17295042316122606504398054547289,
8962 3.2099992681699613513775259670214,
8963 3.24662674946606137764916854570219,
8964 3.28284687953866689817670991319787,
8965 3.31867291347259485044591136879087,
8966 3.35411740487202127264475726990106,
8967 3.38919225660177218727305224515862,
8968 3.42390876691942143189170489271753,
8969 3.45827767149820230182596660024454,
8970 3.49230918177808483937957161007792,
8971 3.5260130200285724149540352829756,
8972 3.55939845146044235497103883695448,
8973 3.59247431368364585025958062194665,
8974 3.62524904377393592090180712976368,
8975 3.65773070318071087226169680450936,
8976 3.68992700068237648299565823810245,
8977 3.72184531357268220291630708234186 };
8979 /* System generated locals */
8980 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
8983 /* Local variables */
8984 static integer ncut, i__;
8985 static doublereal bidon, error;
8986 static integer ia, nd;
8987 static doublereal bid, eps1;
8990 /* ***********************************************************************
8995 /* Lower the degree of a curve defined on (-1,1) in the direction of */
8996 /* Legendre with a given precision. */
9000 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
9002 /* INPUT ARGUMENTS : */
9003 /* ------------------ */
9004 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9005 /* NDIMEN : Dimension of the space. */
9006 /* NCOEFF : Degree +1 of the polynom. */
9007 /* EPSI3D : Precision required for the approximation. */
9008 /* CRVLGD : The curve the degree which of will be lowered. */
9010 /* OUTPUT ARGUMENTS : */
9011 /* ------------------- */
9012 /* YCVMAX : Auxiliary table (error max on each dimension).
9014 /* EPSTRC : Precision of the approximation. */
9015 /* NCFNEW : Degree +1 of the resulting polynom. */
9017 /* COMMONS USED : */
9018 /* ---------------- */
9020 /* REFERENCES CALLED : */
9021 /* ----------------------- */
9023 /* DESCRIPTION/NOTES/LIMITATIONS : */
9024 /* ----------------------------------- */
9026 /* ***********************************************************************
9030 /* Parameter adjustments */
9032 crvlgd_dim1 = *ncofmx;
9033 crvlgd_offset = crvlgd_dim1 + 1;
9034 crvlgd -= crvlgd_offset;
9040 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9044 /* Init for calculation of error. */
9046 for (i__ = 1; i__ <= i__1; ++i__) {
9053 /* Cutting of coefficients. */
9056 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9059 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9060 /* Factor of renormalization. */
9061 bidon = xmaxj[i__ - ncut];
9063 for (nd = 1; nd <= i__2; ++nd) {
9064 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], abs(d__1)) *
9068 /* One stops to cut if the norm becomes too great. */
9069 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9070 if (error > *epsi3d) {
9075 /* --- Max error cumulated when the I-th coeff is removed. */
9082 /* ------- Cutting of zero coeffs of interpolation (RBD) -------
9086 if (*ncfnew == ia) {
9087 AdvApp2Var_MathBase::mmeps1_(&eps1);
9088 for (i__ = ia; i__ >= 2; --i__) {
9091 for (nd = 1; nd <= i__1; ++nd) {
9092 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], abs(d__1));
9101 /* --- If all coeffs can be removed, this is a point. */
9105 /* --------------------------------- End --------------------------------
9112 //=======================================================================
9113 //function : mmtrpj4_
9115 //=======================================================================
9116 int mmtrpj4_(integer *ncofmx,
9125 /* Initialized data */
9127 static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
9128 1.05299572648705464724876659688996,
9129 1.0949715351434178709281698645813,
9130 1.15078388379719068145021100764647,
9131 1.2094863084718701596278219811869,
9132 1.26806623151369531323304177532868,
9133 1.32549784426476978866302826176202,
9134 1.38142537365039019558329304432581,
9135 1.43575531950773585146867625840552,
9136 1.48850442653629641402403231015299,
9137 1.53973611681876234549146350844736,
9138 1.58953193485272191557448229046492,
9139 1.63797820416306624705258190017418,
9140 1.68515974143594899185621942934906,
9141 1.73115699602477936547107755854868,
9142 1.77604489805513552087086912113251,
9143 1.81989256661534438347398400420601,
9144 1.86276344480103110090865609776681,
9145 1.90471563564740808542244678597105,
9146 1.94580231994751044968731427898046,
9147 1.98607219357764450634552790950067,
9148 2.02556989246317857340333585562678,
9149 2.06433638992049685189059517340452,
9150 2.10240936014742726236706004607473,
9151 2.13982350649113222745523925190532,
9152 2.17661085564771614285379929798896,
9153 2.21280102016879766322589373557048,
9154 2.2484214321456956597803794333791,
9155 2.28349755104077956674135810027654,
9156 2.31805304852593774867640120860446,
9157 2.35210997297725685169643559615022,
9158 2.38568889602346315560143377261814,
9159 2.41880904328694215730192284109322,
9160 2.45148841120796359750021227795539,
9161 2.48374387161372199992570528025315,
9162 2.5155912654873773953959098501893,
9163 2.54704548720896557684101746505398,
9164 2.57812056037881628390134077704127,
9165 2.60882970619319538196517982945269,
9166 2.63918540521920497868347679257107,
9167 2.66919945330942891495458446613851,
9168 2.69888301230439621709803756505788,
9169 2.72824665609081486737132853370048,
9170 2.75730041251405791603760003778285,
9171 2.78605380158311346185098508516203,
9172 2.81451587035387403267676338931454,
9173 2.84269522483114290814009184272637,
9174 2.87060005919012917988363332454033,
9175 2.89823818258367657739520912946934,
9176 2.92561704377132528239806135133273,
9177 2.95274375377994262301217318010209,
9178 2.97962510678256471794289060402033,
9179 3.00626759936182712291041810228171,
9180 3.03267744830655121818899164295959,
9181 3.05886060707437081434964933864149 };
9183 /* System generated locals */
9184 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
9187 /* Local variables */
9188 static integer ncut, i__;
9189 static doublereal bidon, error;
9190 static integer ia, nd;
9191 static doublereal bid, eps1;
9195 /* ***********************************************************************
9200 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
9201 /* Legendre with a given precision. */
9205 /* LEGENDRE, POLYGON, TRONCATION, CURVE, SMOOTHING. */
9207 /* INPUT ARGUMENTS : */
9208 /* ------------------ */
9209 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9210 /* NDIMEN : Dimension of the space. */
9211 /* NCOEFF : Degree +1 of the polynom. */
9212 /* EPSI3D : Precision required for the approximation. */
9213 /* CRVLGD : The curve which wishes to lower the degree. */
9215 /* OUTPUT ARGUMENTS : */
9216 /* ------------------- */
9217 /* YCVMAX : Auxiliary table (max error on each dimension).
9219 /* EPSTRC : Precision of the approximation. */
9220 /* NCFNEW : Degree +1 of the resulting polynom. */
9222 /* COMMONS USED : */
9223 /* ---------------- */
9225 /* REFERENCES CALLED : */
9226 /* ----------------------- */
9228 /* DESCRIPTION/NOTES/LIMITATIONS : */
9229 /* ----------------------------------- */
9231 /* ***********************************************************************
9235 /* Parameter adjustments */
9237 crvlgd_dim1 = *ncofmx;
9238 crvlgd_offset = crvlgd_dim1 + 1;
9239 crvlgd -= crvlgd_offset;
9245 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9249 /* Init for error calculation. */
9251 for (i__ = 1; i__ <= i__1; ++i__) {
9258 /* Cutting of coefficients. */
9261 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9264 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9265 /* Factor of renormalization. */
9266 bidon = xmaxj[i__ - ncut];
9268 for (nd = 1; nd <= i__2; ++nd) {
9269 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], abs(d__1)) *
9273 /* Stop cutting if the norm becomes too great. */
9274 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9275 if (error > *epsi3d) {
9280 /* -- Error max cumulated when the I-eme coeff is removed. */
9287 /* ------- Cutting of zero coeffs of the pole of interpolation (RBD) -------
9291 if (*ncfnew == ia) {
9292 AdvApp2Var_MathBase::mmeps1_(&eps1);
9293 for (i__ = ia; i__ >= 2; --i__) {
9296 for (nd = 1; nd <= i__1; ++nd) {
9297 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], abs(d__1));
9306 /* --- If all coeffs can be removed, this is a point. */
9310 /* --------------------------------- End --------------------------------
9317 //=======================================================================
9318 //function : mmtrpj6_
9320 //=======================================================================
9321 int mmtrpj6_(integer *ncofmx,
9331 /* Initialized data */
9333 static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
9334 1.11626917091567929907256116528817,
9335 1.1327140810290884106278510474203,
9336 1.1679452722668028753522098022171,
9337 1.20910611986279066645602153641334,
9338 1.25228283758701572089625983127043,
9339 1.29591971597287895911380446311508,
9340 1.3393138157481884258308028584917,
9341 1.3821288728999671920677617491385,
9342 1.42420414683357356104823573391816,
9343 1.46546895108549501306970087318319,
9344 1.50590085198398789708599726315869,
9345 1.54550385142820987194251585145013,
9346 1.58429644271680300005206185490937,
9347 1.62230484071440103826322971668038,
9348 1.65955905239130512405565733793667,
9349 1.69609056468292429853775667485212,
9350 1.73193098017228915881592458573809,
9351 1.7671112206990325429863426635397,
9352 1.80166107681586964987277458875667,
9353 1.83560897003644959204940535551721,
9354 1.86898184653271388435058371983316,
9355 1.90180515174518670797686768515502,
9356 1.93410285411785808749237200054739,
9357 1.96589749778987993293150856865539,
9358 1.99721027139062501070081653790635,
9359 2.02806108474738744005306947877164,
9360 2.05846864831762572089033752595401,
9361 2.08845055210580131460156962214748,
9362 2.11802334209486194329576724042253,
9363 2.14720259305166593214642386780469,
9364 2.17600297710595096918495785742803,
9365 2.20443832785205516555772788192013,
9366 2.2325216999457379530416998244706,
9367 2.2602654243075083168599953074345,
9368 2.28768115912702794202525264301585,
9369 2.3147799369092684021274946755348,
9370 2.34157220782483457076721300512406,
9371 2.36806787963276257263034969490066,
9372 2.39427635443992520016789041085844,
9373 2.42020656255081863955040620243062,
9374 2.44586699364757383088888037359254,
9375 2.47126572552427660024678584642791,
9376 2.49641045058324178349347438430311,
9377 2.52130850028451113942299097584818,
9378 2.54596686772399937214920135190177,
9379 2.5703922285006754089328998222275,
9380 2.59459096001908861492582631591134,
9381 2.61856915936049852435394597597773,
9382 2.64233265984385295286445444361827,
9383 2.66588704638685848486056711408168,
9384 2.68923766976735295746679957665724,
9385 2.71238965987606292679677228666411 };
9387 /* System generated locals */
9388 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
9391 /* Local variables */
9392 static integer ncut, i__;
9393 static doublereal bidon, error;
9394 static integer ia, nd;
9395 static doublereal bid, eps1;
9399 /* ***********************************************************************
9404 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
9405 /* Legendre to a given precision. */
9409 /* LEGENDRE,POLYGON,TRUNCATION,CURVE,SMOOTHING. */
9411 /* INPUT ARGUMENTS : */
9412 /* ------------------ */
9413 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9414 /* NDIMEN : Dimension of the space. */
9415 /* NCOEFF : Degree +1 of the polynom. */
9416 /* EPSI3D : Precision required for the approximation. */
9417 /* CRVLGD : The curve the degree which of will be lowered. */
9419 /* OUTPUT ARGUMENTS : */
9420 /* ------------------- */
9421 /* YCVMAX : Auxiliary table (max error on each dimension).
9422 /* EPSTRC : Precision of the approximation. */
9423 /* NCFNEW : Degree +1 of the resulting polynom. */
9425 /* COMMONS USED : */
9426 /* ---------------- */
9428 /* REFERENCES CALLED : */
9429 /* ----------------------- */
9431 /* DESCRIPTION/NOTES/LIMITATIONS : */
9432 /* ----------------------------------- */
9434 /* ***********************************************************************
9438 /* Parameter adjustments */
9440 crvlgd_dim1 = *ncofmx;
9441 crvlgd_offset = crvlgd_dim1 + 1;
9442 crvlgd -= crvlgd_offset;
9448 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9452 /* Init for error calculation. */
9454 for (i__ = 1; i__ <= i__1; ++i__) {
9461 /* Cutting of coefficients. */
9464 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9467 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9468 /* Factor of renormalization. */
9469 bidon = xmaxj[i__ - ncut];
9471 for (nd = 1; nd <= i__2; ++nd) {
9472 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], abs(d__1)) *
9476 /* Stop cutting if the norm becomes too great. */
9477 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9478 if (error > *epsi3d) {
9483 /* --- Max error cumulated when the I-th coeff is removed. */
9490 /* ------- Cutting of zero coeff. of the pole of interpolation (RBD) -------
9494 if (*ncfnew == ia) {
9495 AdvApp2Var_MathBase::mmeps1_(&eps1);
9496 for (i__ = ia; i__ >= 2; --i__) {
9499 for (nd = 1; nd <= i__1; ++nd) {
9500 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], abs(d__1));
9509 /* --- If all coeffs can be removed, this is a point. */
9513 /* --------------------------------- End --------------------------------
9520 //=======================================================================
9521 //function : AdvApp2Var_MathBase::mmtrpjj_
9523 //=======================================================================
9524 int AdvApp2Var_MathBase::mmtrpjj_(integer *ncofmx,
9534 /* System generated locals */
9535 integer crvlgd_dim1, crvlgd_offset;
9537 /* Local variables */
9541 /* ***********************************************************************
9546 /* Lower the degree of a curve defined on (-1,1) in the direction of */
9547 /* Legendre with a given precision. */
9551 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
9553 /* INPUT ARGUMENTS : */
9554 /* ------------------ */
9555 /* NCOFMX : Max Nb coeff. of the curve (dimensioning). */
9556 /* NDIMEN : Dimension of the space. */
9557 /* NCOEFF : Degree +1 of the polynom. */
9558 /* EPSI3D : Precision required for the approximation. */
9559 /* IORDRE : Order of continuity at the extremities. */
9560 /* CRVLGD : The curve the degree which of should be lowered. */
9562 /* OUTPUT ARGUMENTS : */
9563 /* ------------------- */
9564 /* ERRMAX : Precision of the approximation. */
9565 /* NCFNEW : Degree +1 of the resulting polynom. */
9567 /* COMMONS USED : */
9568 /* ---------------- */
9570 /* REFERENCES CALLED : */
9571 /* ----------------------- */
9573 /* DESCRIPTION/NOTES/LIMITATIONS : */
9574 /* ----------------------------------- */
9576 /* ***********************************************************************
9580 /* Parameter adjustments */
9582 crvlgd_dim1 = *ncofmx;
9583 crvlgd_offset = crvlgd_dim1 + 1;
9584 crvlgd -= crvlgd_offset;
9587 ia = (*iordre + 1) << 1;
9590 mmtrpj0_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9591 ycvmax[1], errmax, ncfnew);
9592 } else if (ia == 2) {
9593 mmtrpj2_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9594 ycvmax[1], errmax, ncfnew);
9595 } else if (ia == 4) {
9596 mmtrpj4_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9597 ycvmax[1], errmax, ncfnew);
9599 mmtrpj6_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9600 ycvmax[1], errmax, ncfnew);
9603 /* ------------------------ End -----------------------------------------
9609 //=======================================================================
9610 //function : AdvApp2Var_MathBase::mmunivt_
9612 //=======================================================================
9613 int AdvApp2Var_MathBase::mmunivt_(integer *ndimen,
9620 static doublereal c_b2 = 10.;
9622 /* System generated locals */
9626 /* Local variables */
9627 static integer nchif, iunit, izero;
9628 static doublereal vnorm;
9630 static doublereal bid;
9631 static doublereal eps0;
9636 /* ***********************************************************************
9641 /* CALCULATE THE NORMAL VECTOR BASING ON ANY VECTOR */
9642 /* WITH PRECISION GIVEN BY THE USER. */
9646 /* ALL, MATH_ACCES :: */
9647 /* VECTEUR&, NORMALISATION, &VECTEUR */
9649 /* INPUT ARGUMENTS : */
9650 /* ------------------ */
9651 /* NDIMEN : DIMENSION OF THE SPACE */
9652 /* VECTOR : VECTOR TO BE NORMED */
9653 /* EPSILN : EPSILON BELOW WHICH IT IS CONSIDERED THAT THE */
9654 /* NORM OF THE VECTOR IS NULL. IF EPSILN<=0, A DEFAULT VALUE */
9655 /* IS IMPOSED (10.D-17 ON VAX). */
9657 /* OUTPUT ARGUMENTS : */
9658 /* ------------------- */
9659 /* VECNRM : NORMED VECTOR */
9660 /* IERCOD 101 : THE VECTOR IS NULL UP TO EPSILN. */
9663 /* COMMONS USED : */
9664 /* ---------------- */
9666 /* REFERENCES CALLED : */
9667 /* ----------------------- */
9669 /* DESCRIPTION/NOTES/LIMITATIONS : */
9670 /* ----------------------------------- */
9671 /* VECTOR and VECNRM can be identic. */
9673 /* The norm of vector is calculated and each component is divided by
9674 /* this norm. After this it is checked if all componentes of the */
9675 /* vector except for one cost 0 with machine precision. In */
9676 /* this case the quasi-null components are set to 0.D0. */
9678 /* ***********************************************************************
9682 /* Parameter adjustments */
9689 /* -------- Precision by default : zero machine 10.D-17 on Vax ------
9692 AdvApp2Var_SysBase::maovsr8_(&nchif);
9693 if (*epsiln <= 0.) {
9695 eps0 = AdvApp2Var_MathBase::pow__di(&c_b2, &i__1);
9700 /* ------------------------- Calculation of the norm --------------------
9703 vnorm = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vector[1]);
9704 if (vnorm <= eps0) {
9705 AdvApp2Var_SysBase::mvriraz_((integer *)ndimen, (char *)&vecnrm[1]);
9710 /* ---------------------- Calculation of the vector norm ---------------
9714 i__1 = (-nchif - 1) / 2;
9715 eps0 = AdvApp2Var_MathBase::pow__di(&c_b2, &i__1);
9717 for (ii = 1; ii <= i__1; ++ii) {
9718 vecnrm[ii] = vector[ii] / vnorm;
9719 if ((d__1 = vecnrm[ii], abs(d__1)) <= eps0) {
9727 /* ------ Case when all coordinates except for one are almost null ----
9729 /* ------------- then one of coordinates costs 1.D0 or -1.D0 --------
9732 if (izero == *ndimen - 1) {
9733 bid = vecnrm[iunit];
9735 for (ii = 1; ii <= i__1; ++ii) {
9742 vecnrm[iunit] = -1.;
9746 /* -------------------------------- The end -----------------------------
9753 //=======================================================================
9754 //function : AdvApp2Var_MathBase::mmveps3_
9756 //=======================================================================
9757 int AdvApp2Var_MathBase::mmveps3_(doublereal *eps03)
9759 /* Initialized data */
9761 static char nomprg[8+1] = "MMEPS1 ";
9767 /************************************************************************
9772 /* Extraction of EPS1 from COMMON MPRCSN. */
9776 /* MPRCSN,PRECISON,EPS3. */
9778 /* INPUT ARGUMENTS : */
9779 /* ------------------ */
9782 /* OUTPUT ARGUMENTS : */
9783 /* ------------------- */
9784 /* EPS3 : space zero of the denominator (10**-9) */
9785 /* EPS3 should value 10**-15 */
9787 /* COMMONS USED : */
9788 /* ---------------- */
9790 /* REFERENCES CALLED : */
9791 /* ----------------------- */
9793 /* DESCRIPTION/NOTES/LIMITATIONS : */
9794 /* ----------------------------------- */
9797 /* ***********************************************************************
9802 /* ***********************************************************************
9807 /* GIVES TOLERANCES OF NULLITY IN STRIM */
9808 /* AND LIMITS OF ITERATIVE PROCESSES */
9810 /* GENERAL CONTEXT, MODIFIABLE BY THE UTILISER */
9814 /* PARAMETER , TOLERANCE */
9816 /* DESCRIPTION/NOTES/LIMITATIONS : */
9817 /* ----------------------------------- */
9818 /* INITIALISATION : PROFILE , **VIA MPRFTX** AT INPUT IN STRIM*/
9819 /* LOADING OF DEFAULT VALUES OF THE PROFILE IN MPRFTX AT INPUT*/
9820 /* IN STRIM. THEY ARE PRESERVED IN THE LOCAL VARIABLES OF MPRFTX */
9822 /* RESET DEFAULT VALUES : MDFINT */
9823 /* MODIFICATION INTERACTIVE BY THE USER : MDBINT */
9825 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
9826 /* MEPSPB ... EPS3,EPS4 */
9827 /* MEPSLN ... EPS2, NITERM , NITERR */
9828 /* MEPSNR ... EPS2 , NITERM */
9829 /* MITERR ... NITERR */
9832 /* ***********************************************************************
9835 /* NITERM : MAX NB OF ITERATIONS */
9836 /* NITERR : NB OF RAPID ITERATIONS */
9837 /* EPS1 : TOLERANCE OF 3D NULL DISTANCE */
9838 /* EPS2 : TOLERANCE OF ZERO PARAMETRIC DISTANCE */
9839 /* EPS3 : TOLERANCE TO AVOID DIVISION BY 0.. */
9840 /* EPS4 : TOLERANCE ANGULAR */
9844 /* ***********************************************************************
9847 ibb = AdvApp2Var_SysBase::mnfndeb_();
9849 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
9852 *eps03 = mmprcsn_.eps3;
9857 //=======================================================================
9858 //function : AdvApp2Var_MathBase::mmvncol_
9860 //=======================================================================
9861 int AdvApp2Var_MathBase::mmvncol_(integer *ndimen,
9867 /* System generated locals */
9870 /* Local variables */
9871 static logical ldbg;
9873 static doublereal vaux1[3], vaux2[3];
9874 static logical colin;
9875 static doublereal valaux;
9879 /* ***********************************************************************
9884 /* CALCULATE A VECTOR NON-COLINEAR TO A GIVEN NON-NULL VECTOR */
9888 /* PUBLIC, VECTOR, FREE */
9890 /* INPUT ARGUMENTS : */
9891 /* -------------------- */
9892 /* ndimen : dimension of the space */
9893 /* vecin : input vector */
9895 /* OUTPUT ARGUMENTS : */
9896 /* --------------------- */
9898 /* vecout : vector non colinear to vecin */
9900 /* COMMONS USED : */
9901 /* ------------------ */
9904 /* REFERENCES CALLED : */
9905 /* --------------------- */
9908 /* DESCRIPTION/NOTES/LIMITATIONS : */
9909 /* ----------------------------------- */
9911 /* ***********************************************************************
9914 /* ***********************************************************************
9919 /* ***********************************************************************
9921 /* INITIALISATIONS */
9922 /* ***********************************************************************
9925 /* Parameter adjustments */
9930 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
9932 AdvApp2Var_SysBase::mgenmsg_("MMVNCOL", 7L);
9936 /* ***********************************************************************
9939 /* ***********************************************************************
9942 if (*ndimen <= 1 || *ndimen > 3) {
9948 while(d__ <= *ndimen) {
9949 if (vecin[d__] == 0.) {
9954 if (aux == *ndimen) {
9959 for (d__ = 1; d__ <= 3; ++d__) {
9960 vaux1[d__ - 1] = 0.;
9963 for (d__ = 1; d__ <= i__1; ++d__) {
9964 vaux1[d__ - 1] = vecin[d__];
9965 vaux2[d__ - 1] = vecin[d__];
9974 vaux2[d__ - 1] += 1;
9975 valaux = vaux1[1] * vaux2[2] - vaux1[2] * vaux2[1];
9977 valaux = vaux1[2] * vaux2[0] - vaux1[0] * vaux2[2];
9979 valaux = vaux1[0] * vaux2[1] - vaux1[1] * vaux2[0];
9994 for (d__ = 1; d__ <= i__1; ++d__) {
9995 vecout[d__] = vaux2[d__ - 1];
10000 /* ***********************************************************************
10002 /* ERROR PROCESSING */
10003 /* ***********************************************************************
10012 /* ***********************************************************************
10014 /* RETURN CALLING PROGRAM */
10015 /* ***********************************************************************
10021 AdvApp2Var_SysBase::maermsg_("MMVNCOL", iercod, 7L);
10023 AdvApp2Var_SysBase::mgsomsg_("MMVNCOL", 7L);
10028 //=======================================================================
10029 //function : AdvApp2Var_MathBase::mmwprcs_
10031 //=======================================================================
10032 void AdvApp2Var_MathBase::mmwprcs_(doublereal *epsil1,
10033 doublereal *epsil2,
10034 doublereal *epsil3,
10035 doublereal *epsil4,
10042 /* ***********************************************************************
10047 /* ACCESS IN WRITING FOR COMMON MPRCSN */
10053 /* INPUT ARGUMENTS : */
10054 /* -------------------- */
10055 /* EPSIL1 : TOLERANCE OF 3D NULL DISTANCE */
10056 /* EPSIL2 : TOLERANCE OF PARAMETRIC NULL DISTANCE */
10057 /* EPSIL3 : TOLERANCE TO AVOID DIVISION BY 0.. */
10058 /* EPSIL4 : ANGULAR TOLERANCE */
10059 /* NITER1 : MAX NB OF ITERATIONS */
10060 /* NITER2 : NB OF RAPID ITERATIONS */
10062 /* OUTPUT ARGUMENTS : */
10063 /* --------------------- */
10066 /* COMMONS USED : */
10067 /* ------------------ */
10070 /* REFERENCES CALLED : */
10071 /* --------------------- */
10074 /* DESCRIPTION/NOTES/LIMITATIONS : */
10075 /* ----------------------------------- */
10078 /* ***********************************************************************
10081 /* ***********************************************************************
10085 /* ***********************************************************************
10087 /* INITIALIZATIONS */
10088 /* ***********************************************************************
10091 /* ***********************************************************************
10094 /* ***********************************************************************
10097 /* ***********************************************************************
10102 /* GIVES TOLERANCES OF NULLITY IN STRIM */
10103 /* AND LIMITS OF ITERATIVE PROCESSES */
10105 /* GENERAL CONTEXT, MODIFIABLE BY THE UTILISER */
10109 /* PARAMETER , TOLERANCE */
10111 /* DESCRIPTION/NOTES/LIMITATIONS : */
10112 /* ----------------------------------- */
10113 /* INITIALISATION : PROFILE , **VIA MPRFTX** AT INPUT IN STRIM*/
10114 /* LOADING OF DEFAULT VALUES OF THE PROFILE IN MPRFTX AT INPUT*/
10115 /* IN STRIM. THEY ARE PRESERVED IN THE LOCAL VARIABLES OF MPRFTX */
10117 /* RESET DEFAULT VALUES : MDFINT */
10118 /* MODIFICATION INTERACTIVE BY THE USER : MDBINT */
10120 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
10121 /* MEPSPB ... EPS3,EPS4 */
10122 /* MEPSLN ... EPS2, NITERM , NITERR */
10123 /* MEPSNR ... EPS2 , NITERM */
10124 /* MITERR ... NITERR */
10127 /* ***********************************************************************
10130 /* NITERM : MAX NB OF ITERATIONS */
10131 /* NITERR : NB OF RAPID ITERATIONS */
10132 /* EPS1 : TOLERANCE OF 3D NULL DISTANCE */
10133 /* EPS2 : TOLERANCE OF ZERO PARAMETRIC DISTANCE */
10134 /* EPS3 : TOLERANCE TO AVOID DIVISION BY 0.. */
10135 /* EPS4 : TOLERANCE ANGULAR */
10138 /* ***********************************************************************
10140 mmprcsn_.eps1 = *epsil1;
10141 mmprcsn_.eps2 = *epsil2;
10142 mmprcsn_.eps3 = *epsil3;
10143 mmprcsn_.eps4 = *epsil4;
10144 mmprcsn_.niterm = *niter1;
10145 mmprcsn_.niterr = *niter2;
10150 //=======================================================================
10151 //function : AdvApp2Var_MathBase::pow__di
10153 //=======================================================================
10154 doublereal AdvApp2Var_MathBase::pow__di (doublereal *x,
10158 register integer ii ;
10159 doublereal result ;
10162 if ( *n > 0 ) {absolute = *n;}
10163 else {absolute = -*n;}
10164 /* System generated locals */
10165 for(ii = 0 ; ii < absolute ; ii++) {
10169 result = 1.0e0 / result ;
10175 /* **********************************************************************
10180 /* Calculate integer function power not obligatory in the most efficient way ;
10187 /* INPUT ARGUMENTS : */
10188 /* ------------------ */
10189 /* X : argument of X**N */
10192 /* OUTPUT ARGUMENTS : */
10193 /* ------------------- */
10196 /* COMMONS USED : */
10197 /* ---------------- */
10199 /* REFERENCES CALLED : */
10200 /* ----------------------- */
10202 /* DESCRIPTION/NOTES/LIMITATIONS : */
10203 /* ----------------------------------- */
10206 /* ***********************************************************************/
10208 //=======================================================================
10209 //function : pow__ii
10211 //=======================================================================
10212 integer pow__ii(integer *x,
10216 register integer ii ;
10220 if ( *n > 0 ) {absolute = *n;}
10221 else {absolute = -*n;}
10222 /* System generated locals */
10223 for(ii = 0 ; ii < absolute ; ii++) {
10227 result = 1 / result ;
10233 /* **********************************************************************
10235 /* **********************************************************************
10240 /* Calculate integer function power not obligatory in the most efficient way ;
10247 /* INPUT ARGUMENTS : */
10248 /* ------------------ */
10249 /* X : argument of X**N */
10252 /* OUTPUT ARGUMENTS : */
10253 /* ------------------- */
10256 /* COMMONS USED : */
10257 /* ---------------- */
10259 /* REFERENCES CALLED : */
10260 /* ----------------------- */
10262 /* DESCRIPTION/NOTES/LIMITATIONS : */
10263 /* ----------------------------------- */
10266 /* ***********************************************************************/
10268 //=======================================================================
10269 //function : AdvApp2Var_MathBase::msc_
10271 //=======================================================================
10272 doublereal AdvApp2Var_MathBase::msc_(integer *ndimen,
10273 doublereal *vecte1,
10274 doublereal *vecte2)
10277 /* System generated locals */
10279 doublereal ret_val;
10281 /* Local variables */
10282 static integer i__;
10283 static doublereal x;
10287 /************************************************************************
10292 /* Calculate the scalar product of 2 vectors in the space */
10293 /* of dimension NDIMEN. */
10297 /* PRODUCT MSCALAIRE. */
10299 /* INPUT ARGUMENTS : */
10300 /* ------------------ */
10301 /* NDIMEN : Dimension of the space. */
10302 /* VECTE1,VECTE2: Vectors. */
10304 /* OUTPUT ARGUMENTS : */
10305 /* ------------------- */
10307 /* COMMONS USED : */
10308 /* ---------------- */
10310 /* REFERENCES CALLED : */
10311 /* ----------------------- */
10313 /* DESCRIPTION/NOTES/LIMITATIONS : */
10314 /* ----------------------------------- */
10317 /* ***********************************************************************
10321 /* PRODUIT MSCALAIRE */
10322 /* Parameter adjustments */
10326 /* Function Body */
10330 for (i__ = 1; i__ <= i__1; ++i__) {
10331 x += vecte1[i__] * vecte2[i__];
10336 /* ----------------------------------- THE END --------------------------
10342 //=======================================================================
10343 //function : mvcvin2_
10345 //=======================================================================
10346 int mvcvin2_(integer *ncoeff,
10347 doublereal *crvold,
10348 doublereal *crvnew,
10352 /* System generated locals */
10353 integer i__1, i__2;
10355 /* Local variables */
10356 static integer m1jm1, ncfm1, j, k;
10357 static doublereal bid;
10358 static doublereal cij1, cij2;
10362 /************************************************************************
10367 /* INVERSION OF THE PARAMETERS ON CURVE 2D. */
10371 /* CURVE,2D,INVERSION,PARAMETER. */
10373 /* INPUT ARGUMENTS : */
10374 /* ------------------ */
10375 /* NCOEFF : NB OF COEFF OF THE CURVE. */
10376 /* CRVOLD : CURVE OF ORIGIN */
10378 /* OUTPUT ARGUMENTS : */
10379 /* ------------------- */
10380 /* CRVNEW : THE RESULTING CURVE AFTER CHANGE OF T BY 1-T */
10381 /* IERCOD : 0 OK, */
10382 /* 10 NB OF COEFF NULL OR TOO GREAT. */
10384 /* COMMONS USED : */
10385 /* ---------------- */
10388 /* REFERENCES CALLED : */
10389 /* ---------------------- */
10391 /* DESCRIPTION/NOTES/LIMITATIONS : */
10392 /* ----------------------------------- */
10393 /* THE FOLLOWING CALL IS ABSOLUTELY LEGAL : */
10394 /* CALL MVCVIN2(NCOEFF,CURVE,CURVE,IERCOD), THE TABLE CURVE */
10395 /* BECOMES INPUT AND OUTPUT ARGUMENT (RBD). */
10396 /* BECAUSE OF MCCNP, THE NB OF COEFF OF THE CURVE IS LIMITED TO */
10397 /* NDGCNP+1 = 61. */
10400 /* ***********************************************************************
10404 /* **********************************************************************
10409 /* Serves to provide coefficients of the binome (triangle of Pascal). */
10413 /* Coeff of binome from 0 to 60. read only . init par block data */
10415 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
10416 /* ----------------------------------- */
10417 /* The coefficients of the binome form a triangular matrix. */
10418 /* This matrix is completed in table CNP by transposition. */
10419 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
10421 /* Initialization is done by block-data MMLLL09.RES, */
10422 /* created by program MQINICNP.FOR (see the team (AC) ). */
10426 /* **********************************************************************
10431 /* ***********************************************************************
10434 /* Parameter adjustments */
10438 /* Function Body */
10439 if (*ncoeff < 1 || *ncoeff - 1 > 60) {
10446 /* CONSTANT TERM OF THE NEW CURVE */
10451 for (k = 2; k <= i__1; ++k) {
10452 cij1 += crvold[(k << 1) + 1];
10453 cij2 += crvold[(k << 1) + 2];
10457 if (*ncoeff == 1) {
10461 /* INTERMEDIARY POWERS OF THE PARAMETER */
10463 ncfm1 = *ncoeff - 1;
10466 for (j = 2; j <= i__1; ++j) {
10468 cij1 = crvold[(j << 1) + 1];
10469 cij2 = crvold[(j << 1) + 2];
10471 for (k = j + 1; k <= i__2; ++k) {
10472 bid = mmcmcnp_.cnp[k - 1 + (j - 1) * 61];
10473 cij1 += crvold[(k << 1) + 1] * bid;
10474 cij2 += crvold[(k << 1) + 2] * bid;
10476 crvnew[(j << 1) + 1] = cij1 * m1jm1;
10477 crvnew[(j << 1) + 2] = cij2 * m1jm1;
10480 /* TERM OF THE HIGHEST DEGREE */
10482 crvnew[(*ncoeff << 1) + 1] = -crvold[(*ncoeff << 1) + 1] * m1jm1;
10483 crvnew[(*ncoeff << 1) + 2] = -crvold[(*ncoeff << 1) + 2] * m1jm1;
10487 AdvApp2Var_SysBase::maermsg_("MVCVIN2", iercod, 7L);
10492 //=======================================================================
10493 //function : mvcvinv_
10495 //=======================================================================
10496 int mvcvinv_(integer *ncoeff,
10497 doublereal *crvold,
10498 doublereal *crvnew,
10502 /* System generated locals */
10503 integer i__1, i__2;
10505 /* Local variables */
10506 static integer m1jm1, ncfm1, j, k;
10507 static doublereal bid;
10508 //extern /* Subroutine */ int maermsg_();
10509 static doublereal cij1, cij2, cij3;
10512 /* **********************************************************************
10517 /* INVERSION OF THE PARAMETER ON A CURBE 3D (I.E. INVERSION */
10518 /* OF THE DIRECTION OF PARSING). */
10522 /* CURVE,INVERSION,PARAMETER. */
10524 /* INPUT ARGUMENTS : */
10525 /* ------------------ */
10526 /* NCOEFF : NB OF COEFF OF THE CURVE. */
10527 /* CRVOLD : CURVE OF ORIGIN */
10529 /* OUTPUT ARGUMENTS : */
10530 /* ------------------- */
10531 /* CRVNEW : RESULTING CURVE AFTER CHANGE OF T INTO 1-T */
10532 /* IERCOD : 0 OK, */
10533 /* 10 NB OF COEFF NULL OR TOO GREAT. */
10535 /* COMMONS USED : */
10536 /* ---------------- */
10539 /* REFERENCES CALLED : */
10540 /* ---------------------- */
10542 /* DESCRIPTION/NOTES/LIMITATIONS : */
10543 /* ----------------------------------- */
10544 /* THE FOLLOWING CALL IS ABSOLUTELY LEGAL : */
10545 /* CALL MVCVINV(NCOEFF,CURVE,CURVE,IERCOD), TABLE CURVE */
10546 /* BECOMES INPUT AND OUTPUT ARGUMENT (RBD). */
10547 /* THE NUMBER OF COEFF OF THE CURVE IS LIMITED TO NDGCNP+1 = 61 */
10548 /* BECAUSE OF USE OF COMMON MCCNP. */
10550 /* ***********************************************************************
10553 /* **********************************************************************
10558 /* Serves to provide the binomial coefficients (triangle of Pascal). */
10562 /* Binomial Coeff from 0 to 60. read only . init par block data */
10564 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
10565 /* ----------------------------------- */
10566 /* The binomial coefficients form a triangular matrix. */
10567 /* This matrix is completed in table CNP by its transposition. */
10568 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
10570 /* Initialisation is done by block-data MMLLL09.RES, */
10571 /* created by program MQINICNP.FOR (see the team (AC) ). */
10573 /* **********************************************************************
10578 /* ***********************************************************************
10581 /* Parameter adjustments */
10585 /* Function Body */
10586 if (*ncoeff < 1 || *ncoeff - 1 > 60) {
10592 /* CONSTANT TERM OF THE NEW CURVE */
10598 for (k = 2; k <= i__1; ++k) {
10599 cij1 += crvold[k * 3 + 1];
10600 cij2 += crvold[k * 3 + 2];
10601 cij3 += crvold[k * 3 + 3];
10607 if (*ncoeff == 1) {
10611 /* INTERMEDIARY POWER OF THE PARAMETER */
10613 ncfm1 = *ncoeff - 1;
10616 for (j = 2; j <= i__1; ++j) {
10618 cij1 = crvold[j * 3 + 1];
10619 cij2 = crvold[j * 3 + 2];
10620 cij3 = crvold[j * 3 + 3];
10622 for (k = j + 1; k <= i__2; ++k) {
10623 bid = mmcmcnp_.cnp[k - 1 + (j - 1) * 61];
10624 cij1 += crvold[k * 3 + 1] * bid;
10625 cij2 += crvold[k * 3 + 2] * bid;
10626 cij3 += crvold[k * 3 + 3] * bid;
10629 crvnew[j * 3 + 1] = cij1 * m1jm1;
10630 crvnew[j * 3 + 2] = cij2 * m1jm1;
10631 crvnew[j * 3 + 3] = cij3 * m1jm1;
10635 /* TERM OF THE HIGHEST DEGREE */
10637 crvnew[*ncoeff * 3 + 1] = -crvold[*ncoeff * 3 + 1] * m1jm1;
10638 crvnew[*ncoeff * 3 + 2] = -crvold[*ncoeff * 3 + 2] * m1jm1;
10639 crvnew[*ncoeff * 3 + 3] = -crvold[*ncoeff * 3 + 3] * m1jm1;
10642 AdvApp2Var_SysBase::maermsg_("MVCVINV", iercod, 7L);
10646 //=======================================================================
10647 //function : mvgaus0_
10649 //=======================================================================
10650 int mvgaus0_(integer *kindic,
10651 doublereal *urootl,
10652 doublereal *hiltab,
10657 /* System generated locals */
10660 /* Local variables */
10661 static doublereal tamp[40];
10662 static integer ndegl, kg, ii;
10664 /* **********************************************************************
10669 /* Loading of a degree gives roots of LEGENDRE polynom */
10670 /* DEFINED on [-1,1] and weights of Gauss quadrature formulas */
10671 /* (based on corresponding LAGRANGIAN interpolators). */
10672 /* The symmetry relative to 0 is used between [-1,0] and [0,1]. */
10676 /* . VOLUMIC, LEGENDRE, LAGRANGE, GAUSS */
10678 /* INPUT ARGUMENTSE : */
10679 /* ------------------ */
10681 /* KINDIC : Takes values from 1 to 10 depending of the degree */
10682 /* of the used polynom. */
10683 /* The degree of the polynom is equal to 4 k, i.e. 4, 8, */
10684 /* 12, 16, 20, 24, 28, 32, 36 and 40. */
10686 /* OUTPUT ARGUMENTS : */
10687 /* ------------------- */
10689 /* UROOTL : Roots of LEGENDRE polynom in domain [1,0] */
10690 /* given in decreasing order. For domain [-1,0], it is */
10691 /* necessary to take the opposite values. */
10692 /* HILTAB : LAGRANGE interpolators associated to roots. For */
10693 /* opposed roots, interpolatorsare equal. */
10694 /* NBRVAL : Nb of coefficients. Is equal to the half of degree */
10695 /* depending on the symmetry (i.e. 2*KINDIC). */
10697 /* IERCOD : Error code: */
10698 /* < 0 ==> Attention - Warning */
10699 /* =-1 ==> Value of false KINDIC. NBRVAL is forced to 20 */
10701 /* = 0 ==> Everything is OK */
10703 /* COMMON USED : */
10704 /* ---------------- */
10706 /* REFERENCES CALLED : */
10707 /* ------------------- */
10709 /* DESCRIPTION/NOTES/LIMITATIONS : */
10710 /* --------------------------------- */
10711 /* If KINDIC is not correct (i.e < 1 or > 10), the degree is set */
10712 /* to 40 directly (ATTENTION to overload - to avoid it, */
10713 /* preview UROOTL and HILTAB dimensioned at least to 20). */
10715 /* The value of coefficients was calculated with quadruple precision
10716 /* by JJM with help of GD. */
10717 /* Checking of roots was done by GD. */
10719 /* See detailed explications on the listing */
10721 /* **********************************************************************
10725 /* ------------------------------------ */
10726 /* ****** Test validity of KINDIC ** */
10727 /* ------------------------------------ */
10729 /* Parameter adjustments */
10733 /* Function Body */
10736 if (kg < 1 || kg > 10) {
10741 ndegl = *nbrval << 1;
10743 /* ----------------------------------------------------------------------
10745 /* ****** Load NBRVAL positive roots depending on the degree **
10747 /* ----------------------------------------------------------------------
10749 /* ATTENTION : Sign minus (-) in the loop is intentional. */
10751 mmextrl_(&ndegl, tamp);
10753 for (ii = 1; ii <= i__1; ++ii) {
10754 urootl[ii] = -tamp[ii - 1];
10758 /* ------------------------------------------------------------------- */
10759 /* ****** Loading of NBRVAL Gauss weight depending on the degree ** */
10760 /* ------------------------------------------------------------------- */
10762 mmexthi_(&ndegl, tamp);
10764 for (ii = 1; ii <= i__1; ++ii) {
10765 hiltab[ii] = tamp[ii - 1];
10769 /* ------------------------------- */
10770 /* ****** End of sub-program ** */
10771 /* ------------------------------- */
10776 //=======================================================================
10777 //function : mvpscr2_
10779 //=======================================================================
10780 int mvpscr2_(integer *ncoeff,
10781 doublereal *curve2,
10782 doublereal *tparam,
10783 doublereal *pntcrb)
10785 /* System generated locals */
10788 /* Local variables */
10789 static integer ndeg, kk;
10790 static doublereal xxx, yyy;
10794 /* **********************************************************************
10799 /* POSITIONING ON CURVE (NCF,2) IN SPACE OF DIMENSION 2. */
10803 /* TOUS,MATH_ACCES:: COURBE&,POSITIONNEMENT,&POINT. */
10805 /* INPUT ARGUMENTS : */
10806 /* ------------------ */
10807 /* NCOEFF : NUMBER OF COEFFICIENTS OF THE CURVE */
10808 /* CURVE2 : EQUATION OF CURVE 2D */
10809 /* TPARAM : VALUE OF PARAMETER AT GIVEN POINT */
10811 /* OUTPUT ARGUMENTS : */
10812 /* ------------------- */
10813 /* PNTCRB : COORDINATES OF POINT CORRESPONDING TO PARAMETER */
10814 /* TPARAM ON CURVE 2D CURVE2. */
10816 /* COMMONS USED : */
10817 /* ---------------- */
10819 /* REFERENCES CALLED : */
10820 /* ---------------------- */
10822 /* DESCRIPTION/NOTES/LIMITATIONS : */
10823 /* ----------------------------------- */
10824 /* MSCHEMA OF HORNER. */
10827 /* **********************************************************************
10831 /* -------- INITIALIZATIONS AND PROCESSING OF PARTICULAR CASES ----------
10834 /* ---> Cas when NCOEFF > 1 (case STANDARD). */
10835 /* Parameter adjustments */
10839 /* Function Body */
10840 if (*ncoeff >= 2) {
10843 /* ---> Case when NCOEFF <= 1. */
10844 if (*ncoeff <= 0) {
10848 } else if (*ncoeff == 1) {
10849 pntcrb[1] = curve2[3];
10850 pntcrb[2] = curve2[4];
10854 /* -------------------- MSCHEMA OF HORNER (PARTICULAR CASE) --------------
10859 if (*tparam == 1.) {
10863 for (kk = 1; kk <= i__1; ++kk) {
10864 xxx += curve2[(kk << 1) + 1];
10865 yyy += curve2[(kk << 1) + 2];
10869 } else if (*tparam == 0.) {
10870 pntcrb[1] = curve2[3];
10871 pntcrb[2] = curve2[4];
10875 /* ---------------------------- MSCHEMA OF HORNER ------------------------
10877 /* ---> TPARAM is different from 1.D0 and 0.D0. */
10879 ndeg = *ncoeff - 1;
10880 xxx = curve2[(*ncoeff << 1) + 1];
10881 yyy = curve2[(*ncoeff << 1) + 2];
10882 for (kk = ndeg; kk >= 1; --kk) {
10883 xxx = xxx * *tparam + curve2[(kk << 1) + 1];
10884 yyy = yyy * *tparam + curve2[(kk << 1) + 2];
10889 /* ------------------------ RECOVER THE CALCULATED POINT ---------------
10896 /* ------------------------------ THE END -------------------------------
10903 //=======================================================================
10904 //function : mvpscr3_
10906 //=======================================================================
10907 int mvpscr3_(integer *ncoeff,
10908 doublereal *curve3,
10909 doublereal *tparam,
10910 doublereal *pntcrb)
10913 /* System generated locals */
10916 /* Local variables */
10917 static integer ndeg, kk;
10918 static doublereal xxx, yyy, zzz;
10922 /* **********************************************************************
10927 /* POSITIONING ON A CURVE (3,NCF) IN THE SPACE OF DIMENSION 3. */
10931 /* TOUS, MATH_ACCES:: COURBE&,POSITIONNEMENT,&POINT. */
10933 /* INPUT ARGUMENTS : */
10934 /* ------------------ */
10935 /* NCOEFF : NB OF COEFFICIENTS OF THE CURVE */
10936 /* CURVE3 : EQUATION OF CURVE 3D */
10937 /* TPARAM : VALUE OF THE PARAMETER AT THE GIVEN POINT */
10939 /* OUTPUT ARGUMENTS : */
10940 /* ------------------- */
10941 /* PNTCRB : COORDINATES OF THE POINT CORRESPONDING TO PARAMETER */
10942 /* TPARAM ON CURVE 3D CURVE3. */
10944 /* COMMONS USED : */
10945 /* ---------------- */
10947 /* REFERENCES CALLED : */
10948 /* ---------------------- */
10951 /* DESCRIPTION/NOTES/LIMITATIONS : */
10952 /* ----------------------------------- */
10953 /* MSCHEMA OF HORNER. */
10955 /* **********************************************************************
10958 /* **********************************************************************
10962 /* -------- INITIALISATIONS AND PROCESSING OF PARTICULAR CASES ----------
10965 /* ---> Case when NCOEFF > 1 (cas STANDARD). */
10966 /* Parameter adjustments */
10970 /* Function Body */
10971 if (*ncoeff >= 2) {
10974 /* ---> Case when NCOEFF <= 1. */
10975 if (*ncoeff <= 0) {
10980 } else if (*ncoeff == 1) {
10981 pntcrb[1] = curve3[4];
10982 pntcrb[2] = curve3[5];
10983 pntcrb[3] = curve3[6];
10987 /* -------------------- MSCHEMA OF HORNER (PARTICULAR CASE) --------------
10992 if (*tparam == 1.) {
10997 for (kk = 1; kk <= i__1; ++kk) {
10998 xxx += curve3[kk * 3 + 1];
10999 yyy += curve3[kk * 3 + 2];
11000 zzz += curve3[kk * 3 + 3];
11004 } else if (*tparam == 0.) {
11005 pntcrb[1] = curve3[4];
11006 pntcrb[2] = curve3[5];
11007 pntcrb[3] = curve3[6];
11011 /* ---------------------------- MSCHEMA OF HORNER ------------------------
11013 /* ---> Here TPARAM is different from 1.D0 and 0.D0. */
11015 ndeg = *ncoeff - 1;
11016 xxx = curve3[*ncoeff * 3 + 1];
11017 yyy = curve3[*ncoeff * 3 + 2];
11018 zzz = curve3[*ncoeff * 3 + 3];
11019 for (kk = ndeg; kk >= 1; --kk) {
11020 xxx = xxx * *tparam + curve3[kk * 3 + 1];
11021 yyy = yyy * *tparam + curve3[kk * 3 + 2];
11022 zzz = zzz * *tparam + curve3[kk * 3 + 3];
11027 /* ------------------------ RETURN THE CALCULATED POINT ------------------
11035 /* ------------------------------ THE END -------------------------------
11042 //=======================================================================
11043 //function : AdvApp2Var_MathBase::mvsheld_
11045 //=======================================================================
11046 int AdvApp2Var_MathBase::mvsheld_(integer *n,
11052 /* System generated locals */
11053 integer dtab_dim1, dtab_offset, i__1, i__2;
11055 /* Local variables */
11056 static integer incr;
11057 static doublereal dsave;
11058 static integer i3, i4, i5, incrp1;
11061 /************************************************************************
11066 /* PARSING OF COLUMNS OF TABLE OF REAL*8 BY SHELL METHOD*/
11067 /* (IN INCREASING ORDER) */
11071 /* POINT-ENTRY, PARSING, SHELL */
11073 /* INPUT ARGUMENTS : */
11074 /* ------------------ */
11075 /* N : NUMBER OF COLUMNS OF THE TABLE */
11076 /* IS : NUMBER OF LINE OF THE TABLE */
11077 /* DTAB : TABLE OF REAL*8 TO BE PARSED */
11078 /* ICLE : POSITION OF THE KEY ON THE COLUMN */
11080 /* OUTPUT ARGUMENTS : */
11081 /* ------------------- */
11082 /* DTAB : PARSED TABLE */
11084 /* COMMONS USED : */
11085 /* ---------------- */
11088 /* REFERENCES CALLED : */
11089 /* ---------------------- */
11092 /* DESCRIPTION/NOTES/LIMITATIONS : */
11093 /* ----------------------------------- */
11094 /* CLASSIC SHELL METHOD : PARSING BY SERIES */
11095 /* Declaration DTAB(IS, 1) corresponds to DTAB(IS, *) */
11097 /* ***********************************************************************
11101 /* Parameter adjustments */
11103 dtab_offset = dtab_dim1 + 1;
11104 dtab -= dtab_offset;
11106 /* Function Body */
11110 /* ------------------------ */
11112 /* INITIALIZATION OF THE SEQUENCE OF INCREMENTS */
11113 /* FIND THE GREATEST INCREMENT SO THAT INCR < N/9 */
11117 if (incr >= *n / 9) {
11120 /* ----------------------------- */
11121 incr = incr * 3 + 1;
11124 /* LOOP ON INCREMENTS TILL INCR = 1 */
11125 /* PARSING BY SERIES DISTANT FROM INCR */
11129 /* ----------------- */
11131 for (i3 = incrp1; i3 <= i__1; ++i3) {
11132 /* ---------------------- */
11134 /* SET ELEMENT I3 AT ITS PLACE IN THE SERIES */
11141 /* ------------------------- */
11142 if (dtab[*icle + i4 * dtab_dim1] <= dtab[*icle + (i4 + incr) *
11148 for (i5 = 1; i5 <= i__2; ++i5) {
11149 /* ------------------ */
11150 dsave = dtab[i5 + i4 * dtab_dim1];
11151 dtab[i5 + i4 * dtab_dim1] = dtab[i5 + (i4 + incr) * dtab_dim1];
11152 dtab[i5 + (i4 + incr) * dtab_dim1] = dsave;
11163 /* PASSAGE TO THE NEXT INCREMENT */
11174 //=======================================================================
11175 //function : AdvApp2Var_MathBase::mzsnorm_
11177 //=======================================================================
11178 doublereal AdvApp2Var_MathBase::mzsnorm_(integer *ndimen,
11179 doublereal *vecteu)
11182 /* System generated locals */
11184 doublereal ret_val, d__1, d__2;
11186 /* Local variables */
11187 static doublereal xsom;
11188 static integer i__, irmax;
11192 /* ***********************************************************************
11197 /* SERVES to calculate the euclidian norm of a vector : */
11198 /* ____________________________ */
11199 /* Z = V V(1)**2 + V(2)**2 + ... */
11205 /* INPUT ARGUMENTS : */
11206 /* ------------------ */
11207 /* NDIMEN : Dimension of the vector */
11208 /* VECTEU : vector of dimension NDIMEN */
11210 /* OUTPUT ARGUMENTS : */
11211 /* ------------------- */
11212 /* MZSNORM : Value of the euclidian norm of vector VECTEU */
11214 /* COMMONS USED : */
11215 /* ---------------- */
11219 /* REFERENCES CALLED : */
11220 /* ---------------------- */
11222 /* R*8 ABS R*8 SQRT */
11224 /* DESCRIPTION/NOTESS/LIMITATIONS : */
11225 /* ----------------------------------- */
11226 /* To limit the risks of overflow, */
11227 /* the term of the strongest absolute value is factorized : */
11228 /* _______________________ */
11229 /* Z = !V(1)! * V 1 + (V(2)/V(1))**2 + ... */
11232 /* ***********************************************************************
11235 /* ***********************************************************************
11239 /* ***********************************************************************
11242 /* ***********************************************************************
11245 /* ___ Find the strongest absolute value term */
11247 /* Parameter adjustments */
11250 /* Function Body */
11253 for (i__ = 2; i__ <= i__1; ++i__) {
11254 if ((d__1 = vecteu[irmax], abs(d__1)) < (d__2 = vecteu[i__], abs(d__2)
11261 /* ___ Calculate the norme */
11263 if ((d__1 = vecteu[irmax], abs(d__1)) < 1.) {
11266 for (i__ = 1; i__ <= i__1; ++i__) {
11267 /* Computing 2nd power */
11268 d__1 = vecteu[i__];
11269 xsom += d__1 * d__1;
11272 ret_val = sqrt(xsom);
11276 for (i__ = 1; i__ <= i__1; ++i__) {
11277 if (i__ == irmax) {
11280 /* Computing 2nd power */
11281 d__1 = vecteu[i__] / vecteu[irmax];
11282 xsom += d__1 * d__1;
11286 ret_val = (d__1 = vecteu[irmax], abs(d__1)) * sqrt(xsom);
11289 /* ***********************************************************************
11291 /* RETURN CALLING PROGRAM */
11292 /* ***********************************************************************