1 // Copyright (c) 1999-2012 OPEN CASCADE SAS
3 // The content of this file is subject to the Open CASCADE Technology Public
4 // License Version 6.5 (the "License"). You may not use the content of this file
5 // except in compliance with the License. Please obtain a copy of the License
6 // at http://www.opencascade.org and read it completely before using this file.
8 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
9 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
11 // The Original Code and all software distributed under the License is
12 // distributed on an "AS IS" basis, without warranty of any kind, and the
13 // Initial Developer hereby disclaims all such warranties, including without
14 // limitation, any warranties of merchantability, fitness for a particular
15 // purpose or non-infringement. Please see the License for the specific terms
16 // and conditions governing the rights and limitations under the License.
18 // AdvApp2Var_MathBase.cxx
20 #include <AdvApp2Var_SysBase.hxx>
21 #include <AdvApp2Var_Data_f2c.hxx>
22 #include <AdvApp2Var_MathBase.hxx>
23 #include <AdvApp2Var_Data.hxx>
27 int mmchole_(integer *mxcoef,
39 int mmrslss_(integer *mxcoef,
49 int mfac_(doublereal *f,
53 int mmaper0_(integer *ncofmx,
61 int mmaper2_(integer *ncofmx,
70 int mmaper4_(integer *ncofmx,
79 int mmaper6_(integer *ncofmx,
88 int mmarc41_(integer *ndimax,
98 int mmatvec_(integer *nligne,
109 int mmcvstd_(integer *ncofmx,
117 int mmdrvcb_(integer *ideriv,
126 int mmexthi_(integer *ndegre,
130 int mmextrl_(integer *ndegre,
136 int mmherm0_(doublereal *debfin,
140 int mmherm1_(doublereal *debfin,
146 int mmloncv_(integer *ndimax,
155 int mmpojac_(doublereal *tparam,
163 int mmrslw_(integer *normax,
171 int mmtmave_(integer *nligne,
180 int mmtrpj0_(integer *ncofmx,
189 int mmtrpj2_(integer *ncofmx,
199 int mmtrpj4_(integer *ncofmx,
208 int mmtrpj6_(integer *ncofmx,
217 integer pow__ii(integer *x,
221 int mvcvin2_(integer *ncoeff,
227 int mvcvinv_(integer *ncoeff,
233 int mvgaus0_(integer *kindic,
239 int mvpscr2_(integer *ncoeff,
245 int mvpscr3_(integer *ncoeff,
251 doublereal eps1, eps2, eps3, eps4;
252 integer niterm, niterr;
256 doublereal tdebut, tfinal, verifi, cmherm[576];
259 //=======================================================================
260 //function : AdvApp2Var_MathBase::mdsptpt_
262 //=======================================================================
263 int AdvApp2Var_MathBase::mdsptpt_(integer *ndimen,
270 /* System generated locals */
274 /* Local variables */
276 doublereal* differ = 0;
280 /* **********************************************************************
285 /* CALCULATE DISTANCE BETWEEN TWO POINTS */
289 /* DISTANCE,POINT. */
291 /* INPUT ARGUMENTS : */
292 /* ------------------ */
293 /* NDIMEN: Space Dimension. */
294 /* POINT1: Table of coordinates of the 1st point. */
295 /* POINT2: Table of coordinates of the 2nd point. */
297 /* OUTPUT ARGUMENTS : */
298 /* ------------------- */
299 /* DISTAN: Distance between 2 points. */
302 /* ---------------- */
304 /* REFERENCES CALLED : */
305 /* ----------------------- */
307 /* DESCRIPTION/NOTES/LIMITATIONS : */
308 /* ----------------------------------- */
310 /* **********************************************************************
314 /* ***********************************************************************
317 /* ***********************************************************************
320 /* Parameter adjustment */
328 /* ***********************************************************************
331 /* ***********************************************************************
334 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
336 anAdvApp2Var_SysBase.mcrrqst_(&c__8, ndimen, differ, &iofset, &ier);
339 /* --- If allocation is refused, the trivial method is applied. */
345 for (i__ = 1; i__ <= i__1; ++i__) {
346 /* Computing 2nd power */
347 d__1 = point1[i__] - point2[i__];
348 *distan += d__1 * d__1;
350 *distan = sqrt(*distan);
352 /* --- Otherwise MZSNORM is used to minimize the risks of overflow
357 for (i__ = 1; i__ <= i__1; ++i__) {
359 differ[j] = point2[i__] - point1[i__];
362 *distan = AdvApp2Var_MathBase::mzsnorm_(ndimen, &differ[iofset]);
366 /* ***********************************************************************
368 /* RETURN CALLING PROGRAM */
369 /* ***********************************************************************
372 /* --- Dynamic Desallocation */
375 anAdvApp2Var_SysBase.mcrdelt_(&c__8, ndimen, differ, &iofset, &ier);
381 //=======================================================================
384 //=======================================================================
385 int mfac_(doublereal *f,
389 /* System generated locals */
392 /* Local variables */
395 /* FORTRAN CONFORME AU TEXT */
396 /* CALCUL DE MFACTORIEL N */
397 /* Parameter adjustments */
403 for (i__ = 2; i__ <= i__1; ++i__) {
405 f[i__] = i__ * f[i__ - 1];
410 //=======================================================================
411 //function : AdvApp2Var_MathBase::mmapcmp_
413 //=======================================================================
414 int AdvApp2Var_MathBase::mmapcmp_(integer *ndim,
421 /* System generated locals */
422 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
425 /* Local variables */
426 integer ipair, nd, ndegre, impair, ibb, idg;
427 //extern int mgsomsg_();//mgenmsg_(),
429 /* **********************************************************************
434 /* Compression of curve CRVOLD in a table of */
435 /* coeff. of even : CRVNEW(*,0,*) */
436 /* and uneven range : CRVNEW(*,1,*). */
440 /* COMPRESSION,CURVE. */
442 /* INPUT ARGUMENTS : */
443 /* ------------------ */
444 /* NDIM : Space Dimension. */
445 /* NCOFMX : Max nb of coeff. of the curve to compress. */
446 /* NCOEFF : Max nb of coeff. of the compressed curve. */
447 /* CRVOLD : The curve (0:NCOFMX-1,NDIM) to compress. */
449 /* OUTPUT ARGUMENTS : */
450 /* ------------------- */
451 /* CRVNEW : Curve compacted in (0:(NCOEFF-1)/2,0,NDIM) (containing
453 /* even terms) and in (0:(NCOEFF-1)/2,1,NDIM) */
454 /* (containing uneven terms). */
457 /* ---------------- */
459 /* REFERENCES CALLED : */
460 /* ----------------------- */
462 /* DESCRIPTION/NOTES/LIMITATIONS : */
463 /* ----------------------------------- */
464 /* This routine is useful to prepare coefficients of a */
465 /* curve in an orthogonal base (Legendre or Jacobi) before */
466 /* calculating the coefficients in the canonical; base [-1,1] by */
468 /* ***********************************************************************
471 /* Name of the routine */
473 /* Parameter adjustments */
474 crvold_dim1 = *ncofmx;
475 crvold_offset = crvold_dim1;
476 crvold -= crvold_offset;
477 crvnew_dim1 = (*ncoeff - 1) / 2 + 1;
478 crvnew_offset = crvnew_dim1 << 1;
479 crvnew -= crvnew_offset;
482 ibb = AdvApp2Var_SysBase::mnfndeb_();
484 AdvApp2Var_SysBase::mgenmsg_("MMAPCMP", 7L);
487 ndegre = *ncoeff - 1;
489 for (nd = 1; nd <= i__1; ++nd) {
492 for (idg = 0; idg <= i__2; ++idg) {
493 crvnew[idg + (nd << 1) * crvnew_dim1] = crvold[ipair + nd *
502 i__2 = (ndegre - 1) / 2;
503 for (idg = 0; idg <= i__2; ++idg) {
504 crvnew[idg + ((nd << 1) + 1) * crvnew_dim1] = crvold[impair + nd *
515 /* ---------------------------------- The end ---------------------------
519 AdvApp2Var_SysBase::mgsomsg_("MMAPCMP", 7L);
524 //=======================================================================
525 //function : mmaper0_
527 //=======================================================================
528 int mmaper0_(integer *ncofmx,
537 /* System generated locals */
538 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
541 /* Local variables */
546 /* ***********************************************************************
551 /* Calculate the max error of approximation done when */
552 /* only the first NCFNEW coefficients of a curve are preserved.
554 /* Degree NCOEFF-1 written in the base of Legendre (Jacobi */
559 /* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
561 /* INPUT ARGUMENTS : */
562 /* ------------------ */
563 /* NCOFMX : Max. degree of the curve. */
564 /* NDIMEN : Space dimension. */
565 /* NCOEFF : Degree +1 of the curve. */
566 /* CRVLGD : Curve the degree which of should be lowered. */
567 /* NCFNEW : Degree +1 of the resulting polynom. */
569 /* OUTPUT ARGUMENTS : */
570 /* ------------------- */
571 /* YCVMAX : Auxiliary Table (max error on each dimension).
573 /* ERRMAX : Precision of the approximation. */
576 /* ---------------- */
578 /* REFERENCES CALLED : */
579 /* ----------------------- */
581 /* DESCRIPTION/NOTES/LIMITATIONS : */
582 /* ----------------------------------- */
583 /* ***********************************************************************
587 /* ------------------- Init to calculate an error -----------------------
590 /* Parameter adjustments */
592 crvlgd_dim1 = *ncofmx;
593 crvlgd_offset = crvlgd_dim1 + 1;
594 crvlgd -= crvlgd_offset;
598 for (ii = 1; ii <= i__1; ++ii) {
603 /* ------ Minimum that can be reached : Stop at 1 or NCFNEW ------
607 if (*ncfnew + 1 > ncut) {
611 /* -------------- Elimination of high degree coefficients-----------
613 /* ----------- Loop on the series of Legendre: NCUT --> NCOEFF --------
617 for (ii = ncut; ii <= i__1; ++ii) {
618 /* Factor of renormalization (Maximum of Li(t)). */
619 bidon = ((ii - 1) * 2. + 1.) / 2.;
623 for (nd = 1; nd <= i__2; ++nd) {
624 ycvmax[nd] += (d__1 = crvlgd[ii + nd * crvlgd_dim1], advapp_abs(d__1)) *
631 /* -------------- The error is the norm of the vector error ---------------
634 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
636 /* --------------------------------- Fin --------------------------------
642 //=======================================================================
643 //function : mmaper2_
645 //=======================================================================
646 int mmaper2_(integer *ncofmx,
655 /* Initialized data */
657 static doublereal xmaxj[57] = { .9682458365518542212948163499456,
658 .986013297183269340427888048593603,
659 1.07810420343739860362585159028115,
660 1.17325804490920057010925920756025,
661 1.26476561266905634732910520370741,
662 1.35169950227289626684434056681946,
663 1.43424378958284137759129885012494,
664 1.51281316274895465689402798226634,
665 1.5878364329591908800533936587012,
666 1.65970112228228167018443636171226,
667 1.72874345388622461848433443013543,
668 1.7952515611463877544077632304216,
669 1.85947199025328260370244491818047,
670 1.92161634324190018916351663207101,
671 1.98186713586472025397859895825157,
672 2.04038269834980146276967984252188,
673 2.09730119173852573441223706382076,
674 2.15274387655763462685970799663412,
675 2.20681777186342079455059961912859,
676 2.25961782459354604684402726624239,
677 2.31122868752403808176824020121524,
678 2.36172618435386566570998793688131,
679 2.41117852396114589446497298177554,
680 2.45964731268663657873849811095449,
681 2.50718840313973523778244737914028,
682 2.55385260994795361951813645784034,
683 2.59968631659221867834697883938297,
684 2.64473199258285846332860663371298,
685 2.68902863641518586789566216064557,
686 2.73261215675199397407027673053895,
687 2.77551570192374483822124304745691,
688 2.8177699459714315371037628127545,
689 2.85940333797200948896046563785957,
690 2.90044232019793636101516293333324,
691 2.94091151970640874812265419871976,
692 2.98083391718088702956696303389061,
693 3.02023099621926980436221568258656,
694 3.05912287574998661724731962377847,
695 3.09752842783622025614245706196447,
696 3.13546538278134559341444834866301,
697 3.17295042316122606504398054547289,
698 3.2099992681699613513775259670214,
699 3.24662674946606137764916854570219,
700 3.28284687953866689817670991319787,
701 3.31867291347259485044591136879087,
702 3.35411740487202127264475726990106,
703 3.38919225660177218727305224515862,
704 3.42390876691942143189170489271753,
705 3.45827767149820230182596660024454,
706 3.49230918177808483937957161007792,
707 3.5260130200285724149540352829756,
708 3.55939845146044235497103883695448,
709 3.59247431368364585025958062194665,
710 3.62524904377393592090180712976368,
711 3.65773070318071087226169680450936,
712 3.68992700068237648299565823810245,
713 3.72184531357268220291630708234186 };
715 /* System generated locals */
716 integer crvjac_dim1, crvjac_offset, i__1, i__2;
719 /* Local variables */
726 /* ***********************************************************************
731 /* Calculate max approximation error i faite lorsque l' on */
732 /* ne conserve que les premiers NCFNEW coefficients d' une courbe
734 /* de degre NCOEFF-1 ecrite dans la base de Jacobi d' ordre 2. */
738 /* JACOBI, POLYGON, APPROXIMATION, ERROR. */
740 /* INPUT ARGUMENTS : */
741 /* ------------------ */
742 /* NCOFMX : Max. degree of the curve. */
743 /* NDIMEN : Space dimension. */
744 /* NCOEFF : Degree +1 of the curve. */
745 /* CRVLGD : Curve the degree which of should be lowered. */
746 /* NCFNEW : Degree +1 of the resulting polynom. */
748 /* OUTPUT ARGUMENTS : */
749 /* ------------------- */
750 /* YCVMAX : Auxiliary Table (max error on each dimension).
752 /* ERRMAX : Precision of the approximation. */
755 /* ---------------- */
757 /* REFERENCES CALLED : */
758 /* ----------------------- */
759 /* DESCRIPTION/NOTES/LIMITATIONS : */
760 /* ----------------------------------- */
764 /* ------------------ Table of maximums of (1-t2)*Ji(t) ----------------
767 /* Parameter adjustments */
769 crvjac_dim1 = *ncofmx;
770 crvjac_offset = crvjac_dim1 + 1;
771 crvjac -= crvjac_offset;
777 /* ------------------- Init for error calculation -----------------------
781 for (ii = 1; ii <= i__1; ++ii) {
786 /* ------ Min. Degree that can be attained : Stop at 3 or NCFNEW ------
791 i__1 = idec, i__2 = *ncfnew + 1;
792 ncut = advapp_max(i__1,i__2);
794 /* -------------- Removal of coefficients of high degree -----------
796 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
800 for (ii = ncut; ii <= i__1; ++ii) {
801 /* Factor of renormalization. */
802 bidon = xmaxj[ii - idec];
804 for (nd = 1; nd <= i__2; ++nd) {
805 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
812 /* -------------- The error is the norm of the vector error ---------------
815 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
817 /* --------------------------------- Fin --------------------------------
823 /* MAPER4.f -- translated by f2c (version 19960827).
824 You must link the resulting object file with the libraries:
825 -lf2c -lm (in that order)
829 //=======================================================================
830 //function : mmaper4_
832 //=======================================================================
833 int mmaper4_(integer *ncofmx,
841 /* Initialized data */
843 static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
844 1.05299572648705464724876659688996,
845 1.0949715351434178709281698645813,
846 1.15078388379719068145021100764647,
847 1.2094863084718701596278219811869,
848 1.26806623151369531323304177532868,
849 1.32549784426476978866302826176202,
850 1.38142537365039019558329304432581,
851 1.43575531950773585146867625840552,
852 1.48850442653629641402403231015299,
853 1.53973611681876234549146350844736,
854 1.58953193485272191557448229046492,
855 1.63797820416306624705258190017418,
856 1.68515974143594899185621942934906,
857 1.73115699602477936547107755854868,
858 1.77604489805513552087086912113251,
859 1.81989256661534438347398400420601,
860 1.86276344480103110090865609776681,
861 1.90471563564740808542244678597105,
862 1.94580231994751044968731427898046,
863 1.98607219357764450634552790950067,
864 2.02556989246317857340333585562678,
865 2.06433638992049685189059517340452,
866 2.10240936014742726236706004607473,
867 2.13982350649113222745523925190532,
868 2.17661085564771614285379929798896,
869 2.21280102016879766322589373557048,
870 2.2484214321456956597803794333791,
871 2.28349755104077956674135810027654,
872 2.31805304852593774867640120860446,
873 2.35210997297725685169643559615022,
874 2.38568889602346315560143377261814,
875 2.41880904328694215730192284109322,
876 2.45148841120796359750021227795539,
877 2.48374387161372199992570528025315,
878 2.5155912654873773953959098501893,
879 2.54704548720896557684101746505398,
880 2.57812056037881628390134077704127,
881 2.60882970619319538196517982945269,
882 2.63918540521920497868347679257107,
883 2.66919945330942891495458446613851,
884 2.69888301230439621709803756505788,
885 2.72824665609081486737132853370048,
886 2.75730041251405791603760003778285,
887 2.78605380158311346185098508516203,
888 2.81451587035387403267676338931454,
889 2.84269522483114290814009184272637,
890 2.87060005919012917988363332454033,
891 2.89823818258367657739520912946934,
892 2.92561704377132528239806135133273,
893 2.95274375377994262301217318010209,
894 2.97962510678256471794289060402033,
895 3.00626759936182712291041810228171,
896 3.03267744830655121818899164295959,
897 3.05886060707437081434964933864149 };
899 /* System generated locals */
900 integer crvjac_dim1, crvjac_offset, i__1, i__2;
903 /* Local variables */
910 /* ***********************************************************************
915 /* Calculate the max. error of approximation made when */
916 /* only first NCFNEW coefficients of a curve are preserved
918 /* degree NCOEFF-1 is written in the base of Jacobi of order 4. */
921 /* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
923 /* INPUT ARGUMENTS : */
924 /* ------------------ */
925 /* NCOFMX : Max. degree of the curve. */
926 /* NDIMEN : Space dimension. */
927 /* NCOEFF : Degree +1 of the curve. */
928 /* CRVJAC : Curve the degree which of should be lowered. */
929 /* NCFNEW : Degree +1 of the resulting polynom. */
931 /* OUTPUT ARGUMENTS : */
932 /* ------------------- */
933 /* YCVMAX : Auxiliary Table (max error on each dimension).
935 /* ERRMAX : Precision of the approximation. */
938 /* ---------------- */
940 /* REFERENCES CALLED : */
941 /* ----------------------- */
943 /* DESCRIPTION/NOTES/LIMITATIONS : */
946 /* ***********************************************************************
950 /* ---------------- Table of maximums of ((1-t2)2)*Ji(t) ---------------
953 /* Parameter adjustments */
955 crvjac_dim1 = *ncofmx;
956 crvjac_offset = crvjac_dim1 + 1;
957 crvjac -= crvjac_offset;
963 /* ------------------- Init for error calculation -----------------------
967 for (ii = 1; ii <= i__1; ++ii) {
972 /* ------ Min. Degree that can be attained : Stop at 5 or NCFNEW ------
977 i__1 = idec, i__2 = *ncfnew + 1;
978 ncut = advapp_max(i__1,i__2);
980 /* -------------- Removal of high degree coefficients -----------
982 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
986 for (ii = ncut; ii <= i__1; ++ii) {
987 /* Factor of renormalisation. */
988 bidon = xmaxj[ii - idec];
990 for (nd = 1; nd <= i__2; ++nd) {
991 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
998 /* -------------- The error is the norm of the error vector ---------------
1001 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
1003 /* --------------------------------- End --------------------------------
1009 //=======================================================================
1010 //function : mmaper6_
1012 //=======================================================================
1013 int mmaper6_(integer *ncofmx,
1022 /* Initialized data */
1024 static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
1025 1.11626917091567929907256116528817,
1026 1.1327140810290884106278510474203,
1027 1.1679452722668028753522098022171,
1028 1.20910611986279066645602153641334,
1029 1.25228283758701572089625983127043,
1030 1.29591971597287895911380446311508,
1031 1.3393138157481884258308028584917,
1032 1.3821288728999671920677617491385,
1033 1.42420414683357356104823573391816,
1034 1.46546895108549501306970087318319,
1035 1.50590085198398789708599726315869,
1036 1.54550385142820987194251585145013,
1037 1.58429644271680300005206185490937,
1038 1.62230484071440103826322971668038,
1039 1.65955905239130512405565733793667,
1040 1.69609056468292429853775667485212,
1041 1.73193098017228915881592458573809,
1042 1.7671112206990325429863426635397,
1043 1.80166107681586964987277458875667,
1044 1.83560897003644959204940535551721,
1045 1.86898184653271388435058371983316,
1046 1.90180515174518670797686768515502,
1047 1.93410285411785808749237200054739,
1048 1.96589749778987993293150856865539,
1049 1.99721027139062501070081653790635,
1050 2.02806108474738744005306947877164,
1051 2.05846864831762572089033752595401,
1052 2.08845055210580131460156962214748,
1053 2.11802334209486194329576724042253,
1054 2.14720259305166593214642386780469,
1055 2.17600297710595096918495785742803,
1056 2.20443832785205516555772788192013,
1057 2.2325216999457379530416998244706,
1058 2.2602654243075083168599953074345,
1059 2.28768115912702794202525264301585,
1060 2.3147799369092684021274946755348,
1061 2.34157220782483457076721300512406,
1062 2.36806787963276257263034969490066,
1063 2.39427635443992520016789041085844,
1064 2.42020656255081863955040620243062,
1065 2.44586699364757383088888037359254,
1066 2.47126572552427660024678584642791,
1067 2.49641045058324178349347438430311,
1068 2.52130850028451113942299097584818,
1069 2.54596686772399937214920135190177,
1070 2.5703922285006754089328998222275,
1071 2.59459096001908861492582631591134,
1072 2.61856915936049852435394597597773,
1073 2.64233265984385295286445444361827,
1074 2.66588704638685848486056711408168,
1075 2.68923766976735295746679957665724,
1076 2.71238965987606292679677228666411 };
1078 /* System generated locals */
1079 integer crvjac_dim1, crvjac_offset, i__1, i__2;
1082 /* Local variables */
1089 /* ***********************************************************************
1093 /* Calculate the max. error of approximation made when */
1094 /* only first NCFNEW coefficients of a curve are preserved
1096 /* degree NCOEFF-1 is written in the base of Jacobi of order 6. */
1099 /* JACOBI,POLYGON,APPROXIMATION,ERROR. */
1101 /* INPUT ARGUMENTS : */
1102 /* ------------------ */
1103 /* NCOFMX : Max. degree of the curve. */
1104 /* NDIMEN : Space dimension. */
1105 /* NCOEFF : Degree +1 of the curve. */
1106 /* CRVJAC : Curve the degree which of should be lowered. */
1107 /* NCFNEW : Degree +1 of the resulting polynom. */
1109 /* OUTPUT ARGUMENTS : */
1110 /* ------------------- */
1111 /* YCVMAX : Auxiliary Table (max error on each dimension).
1113 /* ERRMAX : Precision of the approximation. */
1115 /* COMMONS USED : */
1116 /* ---------------- */
1118 /* REFERENCES CALLED : */
1119 /* ----------------------- */
1121 /* DESCRIPTION/NOTES/LIMITATIONS : */
1123 /* ***********************************************************************
1127 /* ---------------- Table of maximums of ((1-t2)3)*Ji(t) ---------------
1130 /* Parameter adjustments */
1132 crvjac_dim1 = *ncofmx;
1133 crvjac_offset = crvjac_dim1 + 1;
1134 crvjac -= crvjac_offset;
1140 /* ------------------- Init for error calculation -----------------------
1144 for (ii = 1; ii <= i__1; ++ii) {
1149 /* ------ Min Degree that can be attained : Stop at 3 or NCFNEW ------
1154 i__1 = idec, i__2 = *ncfnew + 1;
1155 ncut = advapp_max(i__1,i__2);
1157 /* -------------- Removal of high degree coefficients -----------
1159 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
1163 for (ii = ncut; ii <= i__1; ++ii) {
1164 /* Factor of renormalization. */
1165 bidon = xmaxj[ii - idec];
1167 for (nd = 1; nd <= i__2; ++nd) {
1168 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
1175 /* -------------- The error is the norm of the vector error ---------------
1178 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
1180 /* --------------------------------- END --------------------------------
1186 //=======================================================================
1187 //function : AdvApp2Var_MathBase::mmaperx_
1189 //=======================================================================
1190 int AdvApp2Var_MathBase::mmaperx_(integer *ncofmx,
1201 /* System generated locals */
1202 integer crvjac_dim1, crvjac_offset;
1204 /* Local variables */
1207 /* **********************************************************************
1211 /* Calculate the max. error of approximation made when */
1212 /* only first NCFNEW coefficients of a curve are preserved
1214 /* degree NCOEFF-1 is written in the base of Jacobi of order IORDRE. */
1217 /* JACOBI,LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
1219 /* INPUT ARGUMENTS : */
1220 /* ------------------ */
1221 /* NCOFMX : Max. degree of the curve. */
1222 /* NDIMEN : Space dimension. */
1223 /* NCOEFF : Degree +1 of the curve. */
1224 /* IORDRE : Order of continuity at the extremities. */
1225 /* CRVJAC : Curve the degree which of should be lowered. */
1226 /* NCFNEW : Degree +1 of the resulting polynom. */
1228 /* OUTPUT ARGUMENTS : */
1229 /* ------------------- */
1230 /* YCVMAX : Auxiliary Table (max error on each dimension).
1232 /* ERRMAX : Precision of the approximation. */
1233 /* IERCOD = 0, OK */
1234 /* = 1, order of constraints (IORDRE) is not within the */
1235 /* autorized values. */
1236 /* COMMONS USED : */
1237 /* ---------------- */
1239 /* REFERENCES CALLED : */
1240 /* ----------------------- */
1242 /* DESCRIPTION/NOTES/LIMITATIONS : */
1243 /* ----------------------------------- */
1244 /* Canceled and replaced MMAPERR. */
1245 /* ***********************************************************************
1249 /* Parameter adjustments */
1251 crvjac_dim1 = *ncofmx;
1252 crvjac_offset = crvjac_dim1 + 1;
1253 crvjac -= crvjac_offset;
1257 /* --> Order of Jacobi polynoms */
1258 jord = ( *iordre + 1) << 1;
1261 mmaper0_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1263 } else if (jord == 2) {
1264 mmaper2_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1266 } else if (jord == 4) {
1267 mmaper4_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1269 } else if (jord == 6) {
1270 mmaper6_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1276 /* ----------------------------------- Fin ------------------------------
1282 //=======================================================================
1283 //function : mmarc41_
1285 //=======================================================================
1286 int mmarc41_(integer *ndimax,
1296 /* System generated locals */
1297 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
1300 /* Local variables */
1302 doublereal tbaux[61];
1308 /* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
1309 /* IMPLICIT INTEGER (I-N) */
1311 /* ***********************************************************************
1316 /* Creation of curve C2(v) defined on (0,1) identic to */
1317 /* curve C1(u) defined on (U0,U1) (change of parameter */
1322 /* LIMITATION, RESTRICTION, CURVE */
1324 /* INPUT ARGUMENTS : */
1325 /* ------------------ */
1326 /* NDIMAX : Space Dimensioning. */
1327 /* NDIMEN : Curve Dimension. */
1328 /* NCOEFF : Nb of coefficients of the curve. */
1329 /* CRVOLD : Curve to be limited. */
1330 /* UPARA0 : Min limit of the interval limiting the curve.
1332 /* UPARA1 : Max limit of the interval limiting the curve.
1335 /* OUTPUT ARGUMENTS : */
1336 /* ------------------- */
1337 /* CRVNEW : Relimited curve, defined on (0,1) and equal to */
1338 /* CRVOLD defined on (U0,U1). */
1339 /* IERCOD : = 0, OK */
1340 /* =10, Nb of coeff. <1 or > 61. */
1342 /* COMMONS USED : */
1343 /* ---------------- */
1344 /* REFERENCES CALLED : */
1345 /* ---------------------- */
1347 /* MAERMSG MCRFILL MVCVIN2 */
1350 /* DESCRIPTION/NOTES/LIMITATIONS : */
1351 /* ----------------------------------- */
1352 /* ---> Algorithm used in this general case is based on the */
1353 /* following principle : */
1354 /* Let S(t) = a0 + a1*t + a2*t**2 + ... of degree NCOEFF-1, and */
1355 /* U(t) = b0 + b1*t, then the coeff. of */
1356 /* S(U(t)) are calculated step by step with help of table TBAUX. */
1357 /* At each step number N (N=2 to NCOEFF), TBAUX(n) contains */
1358 /* the n-th coefficient of U(t)**N for n=1 to N. (RBD) */
1359 /* ---> Reference : KNUTH, 'The Art of Computer Programming', */
1360 /* Vol. 2/'Seminumerical Algorithms', */
1361 /* Ex. 11 p:451 et solution p:562. (RBD) */
1363 /* ---> Removal of the input argument CRVOLD by CRVNEW is */
1364 /* possible, which means that the call : */
1365 /* CALL MMARC41(NDIMAX,NDIMEN,NCOEFF,CURVE,UPARA0,UPARA1 */
1366 /* ,CURVE,IERCOD) */
1367 /* is absolutely LEGAL. (RBD) */
1370 /* **********************************************************************
1373 /* Name of the routine */
1375 /* Auxiliary table of coefficients of (UPARA1-UPARA0)T+UPARA0 */
1376 /* with power N=1 to NCOEFF-1. */
1379 /* Parameter adjustments */
1380 crvnew_dim1 = *ndimax;
1381 crvnew_offset = crvnew_dim1 + 1;
1382 crvnew -= crvnew_offset;
1383 crvold_dim1 = *ndimax;
1384 crvold_offset = crvold_dim1 + 1;
1385 crvold -= crvold_offset;
1389 /* **********************************************************************
1391 /* CASE WHEN PROCESSING CAN'T BE DONE */
1392 /* **********************************************************************
1394 if (*ncoeff > 61 || *ncoeff < 1) {
1398 /* **********************************************************************
1401 /* **********************************************************************
1403 if (*ndimen == *ndimax && *upara0 == 0. && *upara1 == 1.) {
1404 nboct = (*ndimax << 3) * *ncoeff;
1405 AdvApp2Var_SysBase::mcrfill_(&nboct,
1406 &crvold[crvold_offset],
1407 &crvnew[crvnew_offset]);
1410 /* **********************************************************************
1412 /* INVERSION 3D : FAST PROCESSING */
1413 /* **********************************************************************
1415 if (*upara0 == 1. && *upara1 == 0.) {
1416 if (*ndimen == 3 && *ndimax == 3 && *ncoeff <= 21) {
1417 mvcvinv_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
1421 /* ******************************************************************
1423 /* INVERSION 2D : FAST PROCESSING */
1424 /* ******************************************************************
1426 if (*ndimen == 2 && *ndimax == 2 && *ncoeff <= 21) {
1427 mvcvin2_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
1432 /* **********************************************************************
1434 /* GENERAL PROCESSING */
1435 /* **********************************************************************
1437 /* -------------------------- Initializations ---------------------------
1441 for (nd = 1; nd <= i__1; ++nd) {
1442 crvnew[nd + crvnew_dim1] = crvold[nd + crvold_dim1];
1449 tbaux[1] = *upara1 - *upara0;
1451 /* ----------------------- Calculation of coeff. of CRVNEW ------------------
1455 for (ncf = 2; ncf <= i__1; ++ncf) {
1457 /* ------------ Take into account NCF-th coeff. of CRVOLD --------
1461 for (ncj = 1; ncj <= i__2; ++ncj) {
1462 bid = tbaux[ncj - 1];
1464 for (nd = 1; nd <= i__3; ++nd) {
1465 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + ncf *
1472 bid = tbaux[ncf - 1];
1474 for (nd = 1; nd <= i__2; ++nd) {
1475 crvnew[nd + ncf * crvnew_dim1] = crvold[nd + ncf * crvold_dim1] *
1480 /* --------- Calculate (NCF+1) coeff. of ((U1-U0)*t + U0)**(NCF) ---
1483 bid = *upara1 - *upara0;
1484 tbaux[ncf] = tbaux[ncf - 1] * bid;
1485 for (ncj = ncf; ncj >= 2; --ncj) {
1486 tbaux[ncj - 1] = tbaux[ncj - 1] * *upara0 + tbaux[ncj - 2] * bid;
1489 tbaux[0] *= *upara0;
1494 /* -------------- Take into account the last coeff. of CRVOLD -----------
1498 for (ncj = 1; ncj <= i__1; ++ncj) {
1499 bid = tbaux[ncj - 1];
1501 for (nd = 1; nd <= i__2; ++nd) {
1502 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + *ncoeff *
1509 for (nd = 1; nd <= i__1; ++nd) {
1510 crvnew[nd + *ncoeff * crvnew_dim1] = crvold[nd + *ncoeff *
1511 crvold_dim1] * tbaux[*ncoeff - 1];
1515 /* ---------------------------- The end ---------------------------------
1520 AdvApp2Var_SysBase::maermsg_("MMARC41", iercod, 7L);
1526 //=======================================================================
1527 //function : AdvApp2Var_MathBase::mmarcin_
1529 //=======================================================================
1530 int AdvApp2Var_MathBase::mmarcin_(integer *ndimax,
1540 /* System generated locals */
1541 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
1545 /* Local variables */
1548 doublereal tabaux[61];
1557 /* **********************************************************************
1560 /* Creation of curve C2(v) defined on [U0,U1] identic to */
1561 /* curve C1(u) defined on [-1,1] (change of parameter */
1562 /* of a curve) with INVERSION of indices of the resulting table. */
1566 /* GENERALIZED LIMITATION, RESTRICTION, INVERSION, CURVE */
1568 /* INPUT ARGUMENTS : */
1569 /* ------------------ */
1570 /* NDIMAX : Maximum Space Dimensioning. */
1571 /* NDIMEN : Curve Dimension. */
1572 /* NCOEFF : Nb of coefficients of the curve. */
1573 /* CRVOLD : Curve to be limited. */
1574 /* U0 : Min limit of the interval limiting the curve.
1576 /* U1 : Max limit of the interval limiting the curve.
1579 /* OUTPUT ARGUMENTS : */
1580 /* ------------------- */
1581 /* CRVNEW : Relimited curve, defined on [U0,U1] and equal to */
1582 /* CRVOLD defined on [-1,1]. */
1583 /* IERCOD : = 0, OK */
1584 /* =10, Nb of coeff. <1 or > 61. */
1585 /* =13, the requested interval of variation is null. */
1586 /* COMMONS USED : */
1587 /* ---------------- */
1588 /* REFERENCES CALLED : */
1589 /* ---------------------- */
1590 /* DESCRIPTION/NOTES/LIMITATIONS : */
1591 /* ----------------------------------- */
1593 /* **********************************************************************
1596 /* Name of the routine */
1598 /* Auxiliary table of coefficients of X1*T+X0 */
1599 /* with power N=1 to NCOEFF-1. */
1602 /* Parameter adjustments */
1603 crvnew_dim1 = *ndimax;
1604 crvnew_offset = crvnew_dim1 + 1;
1605 crvnew -= crvnew_offset;
1606 crvold_dim1 = *ncoeff;
1607 crvold_offset = crvold_dim1 + 1;
1608 crvold -= crvold_offset;
1611 ibb = AdvApp2Var_SysBase::mnfndeb_();
1613 AdvApp2Var_SysBase::mgenmsg_("MMARCIN", 7L);
1616 /* At zero machine it is tested if the output interval is not null */
1618 AdvApp2Var_MathBase::mmveps3_(&eps3);
1619 if ((d__1 = *u1 - *u0, advapp_abs(d__1)) < eps3) {
1625 /* **********************************************************************
1627 /* CASE WHEN THE PROCESSING IS IMPOSSIBLE */
1628 /* **********************************************************************
1630 if (*ncoeff > 61 || *ncoeff < 1) {
1634 /* **********************************************************************
1636 /* IF NO CHANGE OF THE INTERVAL OF DEFINITION */
1637 /* (ONLY INVERSION OF INDICES OF TABLE CRVOLD) */
1638 /* **********************************************************************
1640 if (*ndim == *ndimax && *u0 == -1. && *u1 == 1.) {
1641 AdvApp2Var_MathBase::mmcvinv_(ndim, ncoeff, ndim, &crvold[crvold_offset], &crvnew[
1645 /* **********************************************************************
1647 /* CASE WHEN THE NEW INTERVAL OF DEFINITION IS [0,1] */
1648 /* **********************************************************************
1650 if (*u0 == 0. && *u1 == 1.) {
1651 mmcvstd_(ncoeff, ndimax, ncoeff, ndim, &crvold[crvold_offset], &
1652 crvnew[crvnew_offset]);
1655 /* **********************************************************************
1657 /* GENERAL PROCESSING */
1658 /* **********************************************************************
1660 /* -------------------------- Initialization ---------------------------
1663 x0 = -(*u1 + *u0) / (*u1 - *u0);
1664 x1 = 2. / (*u1 - *u0);
1666 for (nd = 1; nd <= i__1; ++nd) {
1667 crvnew[nd + crvnew_dim1] = crvold[nd * crvold_dim1 + 1];
1676 /* ----------------------- Calculation of coeff. of CRVNEW ------------------
1680 for (ncf = 2; ncf <= i__1; ++ncf) {
1682 /* ------------ Take into account the NCF-th coeff. of CRVOLD --------
1686 for (ncj = 1; ncj <= i__2; ++ncj) {
1687 bid = tabaux[ncj - 1];
1689 for (nd = 1; nd <= i__3; ++nd) {
1690 crvnew[nd + ncj * crvnew_dim1] += crvold[ncf + nd *
1697 bid = tabaux[ncf - 1];
1699 for (nd = 1; nd <= i__2; ++nd) {
1700 crvnew[nd + ncf * crvnew_dim1] = crvold[ncf + nd * crvold_dim1] *
1705 /* --------- Calculation of (NCF+1) coeff. of [X1*t + X0]**(NCF) --------
1708 tabaux[ncf] = tabaux[ncf - 1] * x1;
1709 for (ncj = ncf; ncj >= 2; --ncj) {
1710 tabaux[ncj - 1] = tabaux[ncj - 1] * x0 + tabaux[ncj - 2] * x1;
1718 /* -------------- Take into account the last coeff. of CRVOLD -----------
1722 for (ncj = 1; ncj <= i__1; ++ncj) {
1723 bid = tabaux[ncj - 1];
1725 for (nd = 1; nd <= i__2; ++nd) {
1726 crvnew[nd + ncj * crvnew_dim1] += crvold[*ncoeff + nd *
1733 for (nd = 1; nd <= i__1; ++nd) {
1734 crvnew[nd + *ncoeff * crvnew_dim1] = crvold[*ncoeff + nd *
1735 crvold_dim1] * tabaux[*ncoeff - 1];
1739 /* ---------------------------- The end ---------------------------------
1744 AdvApp2Var_SysBase::maermsg_("MMARCIN", iercod, 7L);
1747 AdvApp2Var_SysBase::mgsomsg_("MMARCIN", 7L);
1752 //=======================================================================
1753 //function : mmatvec_
1755 //=======================================================================
1756 int mmatvec_(integer *nligne,
1767 /* System generated locals */
1770 /* Local variables */
1772 integer jmin, jmax, i__, j, k;
1777 /* ***********************************************************************
1782 /* Produce vector matrix in form of profile */
1787 /* RESERVE, MATRIX, PRODUCT, VECTOR, PROFILE */
1789 /* INPUT ARGUMENTS : */
1790 /* -------------------- */
1791 /* NLIGNE : Line number of the matrix of constraints */
1792 /* NCOLON : Number of column of the matrix of constraints */
1793 /* GNSTOC: Number of coefficients in the profile of matrix GMATRI */
1795 /* GPOSIT: Table of positioning of terms of storage */
1796 /* GPOSIT(1,I) contains the number of terms-1 on the line I
1797 /* in the profile of the matrix. */
1798 /* GPOSIT(2,I) contains the index of storage of diagonal term*/
1800 /* GPOSIT(3,I) contains the index of column of the first term of */
1801 /* profile of line I */
1802 /* GNSTOC: Number of coefficients in the profile of matrix */
1804 /* GMATRI : Matrix of constraints in form of profile */
1805 /* VECIN : Input vector */
1806 /* DEBLIG : Line indexusing which the vector matrix is calculated */
1808 /* OUTPUT ARGUMENTS */
1809 /* --------------------- */
1810 /* VECOUT : VECTOR PRODUCT */
1812 /* IERCOD : ERROR CODE */
1815 /* COMMONS USED : */
1816 /* ------------------ */
1819 /* REFERENCES CALLED : */
1820 /* --------------------- */
1823 /* DESCRIPTION/NOTES/LIMITATIONS : */
1824 /* ----------------------------------- */
1826 /* ***********************************************************************
1829 /* ***********************************************************************
1834 /* ***********************************************************************
1836 /* INITIALISATIONS */
1837 /* ***********************************************************************
1840 /* Parameter adjustments */
1847 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1849 AdvApp2Var_SysBase::mgenmsg_("MMATVEC", 7L);
1853 /* ***********************************************************************
1856 /* ***********************************************************************
1858 AdvApp2Var_SysBase::mvriraz_(nligne,
1861 for (i__ = *deblig; i__ <= i__1; ++i__) {
1863 jmin = gposit[i__ * 3 + 3];
1864 jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
1865 aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
1867 for (j = jmin; j <= i__2; ++j) {
1869 somme += gmatri[k] * vecin[j];
1871 vecout[i__] = somme;
1880 /* ***********************************************************************
1882 /* ERROR PROCESSING */
1883 /* ***********************************************************************
1889 /* ***********************************************************************
1891 /* RETURN CALLING PROGRAM */
1892 /* ***********************************************************************
1897 /* ___ DESALLOCATION, ... */
1899 AdvApp2Var_SysBase::maermsg_("MMATVEC", iercod, 7L);
1901 AdvApp2Var_SysBase::mgsomsg_("MMATVEC", 7L);
1907 //=======================================================================
1908 //function : mmbulld_
1910 //=======================================================================
1911 int AdvApp2Var_MathBase::mmbulld_(integer *nbcoln,
1917 /* System generated locals */
1918 integer dtabtr_dim1, dtabtr_offset, i__1, i__2;
1920 /* Local variables */
1923 integer nite1, nite2, nchan, i1, i2;
1925 /* ***********************************************************************
1930 /* Parsing of columns of a table of integers in increasing order */
1933 /* POINT-ENTRY, PARSING */
1934 /* INPUT ARGUMENTS : */
1935 /* -------------------- */
1936 /* - NBCOLN : Number of columns in the table */
1937 /* - NBLIGN : Number of lines in the table */
1938 /* - DTABTR : Table of integers to be parsed */
1939 /* - NUMCLE : Position of the key on the column */
1941 /* OUTPUT ARGUMENTS : */
1942 /* --------------------- */
1943 /* - DTABTR : Parsed table */
1945 /* COMMONS USED : */
1946 /* ------------------ */
1949 /* REFERENCES CALLED : */
1950 /* --------------------- */
1953 /* DESCRIPTION/NOTES/LIMITATIONS : */
1954 /* ----------------------------------- */
1955 /* Particularly performant if the table is almost parsed */
1956 /* In the opposite case it is better to use MVSHELD */
1957 /* ***********************************************************************
1960 /* Parameter adjustments */
1961 dtabtr_dim1 = *nblign;
1962 dtabtr_offset = dtabtr_dim1 + 1;
1963 dtabtr -= dtabtr_offset;
1966 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1968 AdvApp2Var_SysBase::mgenmsg_("MMBULLD", 7L);
1974 /* ***********************************************************************
1977 /* ***********************************************************************
1980 /* ---->ALGORITHM in N^2 / 2 additional iteration */
1984 /* ----> Parsing from left to the right */
1988 for (i1 = nite2; i1 <= i__1; ++i1) {
1989 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1 - 1)
1992 for (i2 = 1; i2 <= i__2; ++i2) {
1993 daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
1994 dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
1996 dtabtr[i2 + i1 * dtabtr_dim1] = daux;
2005 /* ----> Parsing from right to the left */
2010 for (i1 = nite1; i1 >= i__1; --i1) {
2011 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1
2012 - 1) * dtabtr_dim1]) {
2014 for (i2 = 1; i2 <= i__2; ++i2) {
2015 daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
2016 dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
2018 dtabtr[i2 + i1 * dtabtr_dim1] = daux;
2032 /* ***********************************************************************
2034 /* ERROR PROCESSING */
2035 /* ***********************************************************************
2038 /* ----> No errors at calling functions, only tests and loops. */
2040 /* ***********************************************************************
2042 /* RETURN CALLING PROGRAM */
2043 /* ***********************************************************************
2049 AdvApp2Var_SysBase::mgsomsg_("MMBULLD", 7L);
2056 //=======================================================================
2057 //function : AdvApp2Var_MathBase::mmcdriv_
2059 //=======================================================================
2060 int AdvApp2Var_MathBase::mmcdriv_(integer *ndimen,
2069 /* System generated locals */
2070 integer courbe_dim1, courbe_offset, crvdrv_dim1, crvdrv_offset, i__1,
2073 /* Local variables */
2075 doublereal mfactk, bid;
2078 /* ***********************************************************************
2083 /* Calculate matrix of a derivate curve of order IDERIV. */
2084 /* with input parameters other than output parameters. */
2089 /* COEFFICIENTS,CURVE,DERIVATE I-EME. */
2091 /* INPUT ARGUMENTS : */
2092 /* ------------------ */
2093 /* NDIMEN : Space dimension (2 or 3 in general) */
2094 /* NCOEFF : Degree +1 of the curve. */
2095 /* COURBE : Table of coefficients of the curve. */
2096 /* IDERIV : Required order of derivation : 1=1st derivate, etc... */
2098 /* OUTPUT ARGUMENTS : */
2099 /* ------------------- */
2100 /* NCOFDV : Degree +1 of the derivative of order IDERIV of the curve. */
2101 /* CRVDRV : Table of coefficients of the derivative of order IDERIV */
2104 /* COMMONS USED : */
2105 /* ---------------- */
2107 /* REFERENCES CALLED : */
2108 /* ----------------------- */
2110 /* DESCRIPTION/NOTES/LIMITATIONS : */
2111 /* ----------------------------------- */
2113 /* ---> It is possible to take as output argument the curve */
2114 /* and the number of coeff passed at input by making : */
2115 /* CALL MMCDRIV(NDIMEN,NCOEFF,COURBE,IDERIV,NCOEFF,COURBE). */
2116 /* After this call, NCOEFF does the number of coeff of the derived */
2117 /* curve the coefficients which of are stored in CURVE. */
2118 /* Attention to the coefficients of CURVE of rank superior to */
2119 /* NCOEFF : they are not set to zero. */
2121 /* ---> Algorithm : */
2122 /* The code below was written basing on the following algorithm:
2125 /* Let P(t) = a1 + a2*t + ... an*t**n. Derivate of order k of P */
2126 /* (containing n-k coefficients) is calculated as follows : */
2128 /* Pk(t) = a(k+1)*CNP(k,k)*k! */
2129 /* + a(k+2)*CNP(k+1,k)*k! * t */
2133 /* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
2134 /* ***********************************************************************
2138 /* -------------- Case when the order of derivative is -------------------
2140 /* ---------------- greater than the degree of the curve ---------------------
2143 /* **********************************************************************
2148 /* Serves to provide the coefficients of binome (Pascal's triangle). */
2152 /* Binomial coeff from 0 to 60. read only . init par block data */
2154 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
2155 /* ----------------------------------- */
2156 /* Binomial coefficients form a triangular matrix. */
2157 /* This matrix is completed in table CNP by its transposition. */
2158 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
2160 /* Initialization is done by block-data MMLLL09.RES, */
2161 /* created by program MQINICNP.FOR). */
2162 /* **********************************************************************
2167 /* ***********************************************************************
2170 /* Parameter adjustments */
2171 crvdrv_dim1 = *ndimen;
2172 crvdrv_offset = crvdrv_dim1 + 1;
2173 crvdrv -= crvdrv_offset;
2174 courbe_dim1 = *ndimen;
2175 courbe_offset = courbe_dim1 + 1;
2176 courbe -= courbe_offset;
2179 if (*ideriv >= *ncoeff) {
2181 for (i__ = 1; i__ <= i__1; ++i__) {
2182 crvdrv[i__ + crvdrv_dim1] = 0.;
2188 /* **********************************************************************
2190 /* General processing */
2191 /* **********************************************************************
2193 /* --------------------- Calculation of Factorial(IDERIV) ------------------
2199 for (i__ = 2; i__ <= i__1; ++i__) {
2204 /* ------------ Calculation of coeff of the derived of order IDERIV ----------
2206 /* ---> Attention : coefficient binomial C(n,m) is represented in */
2207 /* MCCNP by CNP(N+1,M+1). */
2210 for (j = k + 1; j <= i__1; ++j) {
2211 bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
2213 for (i__ = 1; i__ <= i__2; ++i__) {
2214 crvdrv[i__ + (j - k) * crvdrv_dim1] = bid * courbe[i__ + j *
2221 *ncofdv = *ncoeff - *ideriv;
2223 /* -------------------------------- The end -----------------------------
2230 //=======================================================================
2231 //function : AdvApp2Var_MathBase::mmcglc1_
2233 //=======================================================================
2234 int AdvApp2Var_MathBase::mmcglc1_(integer *ndimax,
2247 /* System generated locals */
2248 integer courbe_dim1, courbe_offset, i__1;
2251 /* Local variables */
2253 doublereal tdeb, tfin;
2259 doublereal dif, pas;
2263 /* ***********************************************************************
2268 /* Allows calculating the length of an arc of curve POLYNOMIAL */
2269 /* on an interval [A,B]. */
2273 /* LENGTH,CURVE,GAUSS,PRIVATE. */
2275 /* INPUT ARGUMENTS : */
2276 /* ------------------ */
2277 /* NDIMAX : Max. number of lines of tables */
2278 /* (i.e. max. nb of polynoms). */
2279 /* NDIMEN : Dimension of the space (nb of polynoms). */
2280 /* NCOEFF : Nb of coefficients of the polynom. This is degree + 1.
2282 /* COURBE(NDIMAX,NCOEFF) : Coefficients of the curve. */
2283 /* TDEBUT : Lower limit of the interval of integration for */
2284 /* length calculation. */
2285 /* TFINAL : Upper limit of the interval of integration for */
2286 /* length calculation. */
2287 /* EPSILN : REQIRED precision for length calculation. */
2289 /* OUTPUT ARGUMENTS : */
2290 /* ------------------- */
2291 /* XLONGC : Length of the arc of curve */
2292 /* ERREUR : Precision OBTAINED for the length calculation. */
2293 /* IERCOD : Error code, 0 OK, >0 Serious error. */
2294 /* = 1 Too much iterations, the best calculated resultat */
2295 /* (is almost ERROR) */
2296 /* = 2 Pb MMLONCV (no result) */
2297 /* = 3 NDIM or NCOEFF invalid (no result) */
2299 /* COMMONS USED : */
2300 /* ---------------- */
2302 /* REFERENCES CALLED : */
2303 /* ----------------------- */
2305 /* DESCRIPTION/NOTES/LIMITATIONS : */
2306 /* ----------------------------------- */
2307 /* The polynom is actually a set of polynoms with */
2308 /* coefficients arranged in a table of 2 indices, */
2309 /* each line relative to the polynom. */
2310 /* The polynom is defined by these coefficients ordered */
2311 /* by increasing power of the variable. */
2312 /* All polynoms have the same number of coefficients (the */
2315 /* This program cancels and replaces LENGCV, MLONGC and MLENCV. */
2317 /* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
2320 /* ***********************************************************************
2323 /* Name of the routine */
2326 /* ------------------------ General Initialization ---------------------
2329 /* Parameter adjustments */
2330 courbe_dim1 = *ndimax;
2331 courbe_offset = courbe_dim1 + 1;
2332 courbe -= courbe_offset;
2335 ibb = AdvApp2Var_SysBase::mnfndeb_();
2337 AdvApp2Var_SysBase::mgenmsg_("MMCGLC1", 7L);
2344 /* ------ Test of equity of limits */
2346 if (*tdebut == *tfinal) {
2351 /* ------ Test of the dimension and the number of coefficients */
2353 if (*ndimen <= 0 || *ncoeff <= 0) {
2357 /* ----- Nb of current cutting, nb of iteration, */
2358 /* max nb of iterations */
2365 /* ------ Variation of the nb of intervals */
2366 /* Multiplied by 2 at each iteration */
2369 pas = (*tfinal - *tdebut) / ndec;
2372 /* ------ Loop on all current NDEC intervals */
2375 for (kk = 1; kk <= i__1; ++kk) {
2377 /* ------ Limits of the current integration interval */
2379 tdeb = *tdebut + (kk - 1) * pas;
2381 mmloncv_(ndimax, ndimen, ncoeff, &courbe[courbe_offset], &tdeb, &tfin,
2393 /* ----------------- Test of the maximum number of iterations ------------
2396 /* Test if passes at least once ** */
2405 /* ------ Take into account DIF - Test of convergence */
2408 dif = (d__1 = sottc - oldso, advapp_abs(d__1));
2410 /* ------ If DIF is OK, leave..., otherwise: */
2412 if (dif > *epsiln) {
2414 /* ------ If nb iteration exceeded, leave */
2421 /* ------ Otherwise continue by cutting the initial interval.
2431 /* ------------------------------ THE END -------------------------------
2439 /* ---> PB in MMLONCV */
2445 /* ---> NCOEFF or NDIM invalid. */
2453 AdvApp2Var_SysBase::maermsg_("MMCGLC1", iercod, 7L);
2456 AdvApp2Var_SysBase::mgsomsg_("MMCGLC1", 7L);
2461 //=======================================================================
2462 //function : mmchole_
2464 //=======================================================================
2465 int mmchole_(integer *,//mxcoef,
2474 /* System generated locals */
2475 integer i__1, i__2, i__3;
2478 /* Builtin functions */
2481 /* Local variables */
2483 integer kmin, i__, j, k;
2485 integer ptini, ptcou;
2488 /* ***********************************************************************
2493 /* Produce decomposition of choleski of matrix A in S.S */
2494 /* Calculate inferior triangular matrix S. */
2498 /* RESOLUTION, MFACTORISATION, MATRIX_PROFILE, CHOLESKI */
2500 /* INPUT ARGUMENTS : */
2501 /* -------------------- */
2502 /* MXCOEF : Max number of terms in the hessian profile */
2503 /* DIMENS : Dimension of the problem */
2504 /* AMATRI(MXCOEF) : Coefficients of the matrix profile */
2505 /* APOSIT(1,*) : Distance diagonal-left extremity of the line
2507 /* APOSIT(2,*) : Position of diagonal terms in HESSIE */
2508 /* POSUIV(MXCOEF) : first line inferior not out of profile */
2510 /* OUTPUT ARGUMENTS : */
2511 /* --------------------- */
2512 /* CHOMAT(MXCOEF) : Inferior triangular matrix preserving the */
2513 /* profile of AMATRI. */
2514 /* IERCOD : error code */
2516 /* = 1 : non-defined positive matrix */
2518 /* COMMONS USED : */
2519 /* ------------------ */
2523 /* REFERENCES CALLED : */
2524 /* ---------------------- */
2526 /* DESCRIPTION/NOTES/LIMITATIONS : */
2527 /* ----------------------------------- */
2528 /* DEBUG LEVEL = 4 */
2529 /* ***********************************************************************
2532 /* ***********************************************************************
2537 /* ***********************************************************************
2539 /* INITIALISATIONS */
2540 /* ***********************************************************************
2543 /* Parameter adjustments */
2550 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
2552 AdvApp2Var_SysBase::mgenmsg_("MMCHOLE", 7L);
2556 /* ***********************************************************************
2559 /* ***********************************************************************
2563 for (j = 1; j <= i__1; ++j) {
2565 ptini = aposit[(j << 1) + 2];
2569 for (k = ptini - aposit[(j << 1) + 1]; k <= i__2; ++k) {
2570 /* Computing 2nd power */
2572 somme += d__1 * d__1;
2575 if (amatri[ptini] - somme < 1e-32) {
2578 chomat[ptini] = sqrt(amatri[ptini] - somme);
2582 while(posuiv[ptcou] > 0) {
2584 i__ = posuiv[ptcou];
2585 ptcou = aposit[(i__ << 1) + 2] - (i__ - j);
2587 /* Calculate the sum of S .S for k =1 a j-1 */
2591 i__2 = i__ - aposit[(i__ << 1) + 1], i__3 = j - aposit[(j << 1) +
2593 kmin = advapp_max(i__2,i__3);
2595 for (k = kmin; k <= i__2; ++k) {
2596 somme += chomat[aposit[(i__ << 1) + 2] - (i__ - k)] * chomat[
2597 aposit[(j << 1) + 2] - (j - k)];
2600 chomat[ptcou] = (amatri[ptcou] - somme) / chomat[ptini];
2606 /* ***********************************************************************
2608 /* ERROR PROCESSING */
2609 /* ***********************************************************************
2616 /* ***********************************************************************
2618 /* RETURN CALLING PROGRAM */
2619 /* ***********************************************************************
2624 AdvApp2Var_SysBase::maermsg_("MMCHOLE", iercod, 7L);
2626 AdvApp2Var_SysBase::mgsomsg_("MMCHOLE", 7L);
2632 //=======================================================================
2633 //function : AdvApp2Var_MathBase::mmcvctx_
2635 //=======================================================================
2636 int AdvApp2Var_MathBase::mmcvctx_(integer *ndimen,
2646 /* System generated locals */
2647 integer ctrtes_dim1, ctrtes_offset, crvres_dim1, crvres_offset,
2648 xmatri_dim1, xmatri_offset, tabaux_dim1, tabaux_offset, i__1,
2651 /* Local variables */
2652 integer moup1, nordr;
2654 integer ibb, ncf, ndv;
2658 /* ***********************************************************************
2663 /* Calculate a polynomial curve checking the */
2664 /* passage constraints (interpolation) */
2665 /* from first derivatives, etc... to extremities. */
2666 /* Parameters at the extremities are supposed to be -1 and 1. */
2670 /* ALL, AB_SPECIFI::CONSTRAINTS&,INTERPOLATION,&CURVE */
2672 /* INPUT ARGUMENTS : */
2673 /* ------------------ */
2674 /* NDIMEN : Space Dimension. */
2675 /* NCOFMX : Nb of coeff. of curve CRVRES on each */
2677 /* NDERIV : Order of constraint with derivatives : */
2678 /* 0 --> interpolation simple. */
2679 /* 1 --> interpolation+constraints with 1st. */
2680 /* 2 --> cas (0)+ (1) + " " 2nd derivatives. */
2682 /* CTRTES : Table of constraints. */
2683 /* CTRTES(*,1,*) = contraints at -1. */
2684 /* CTRTES(*,2,*) = contraints at 1. */
2686 /* OUTPUT ARGUMENTS : */
2687 /* ------------------- */
2688 /* CRVRES : Resulting curve defined on (-1,1). */
2689 /* TABAUX : Auxilliary matrix. */
2690 /* XMATRI : Auxilliary matrix. */
2692 /* COMMONS UTILISES : */
2693 /* ---------------- */
2697 /* REFERENCES CALLED : */
2698 /* ---------------------- */
2700 /* MAERMSG R*8 DFLOAT MGENMSG */
2701 /* MGSOMSG MMEPS1 MMRSLW */
2704 /* DESCRIPTION/NOTES/LIMITATIONS : */
2705 /* ----------------------------------- */
2706 /* The polynom (or the curve) is calculated by solving a */
2707 /* system of linear equations. If the imposed degree is great */
2708 /* it is preferable to call a routine based on */
2709 /* Lagrange or Hermite interpolation depending on the case. */
2710 /* (for a high degree the matrix of the system can be badly */
2711 /* conditionned). */
2712 /* This routine returns a curve defined in (-1,1). */
2713 /* In general case, it is necessary to use MCVCTG. */
2715 /* ***********************************************************************
2718 /* Name of the routine */
2721 /* Parameter adjustments */
2722 crvres_dim1 = *ncofmx;
2723 crvres_offset = crvres_dim1 + 1;
2724 crvres -= crvres_offset;
2725 xmatri_dim1 = *nderiv + 1;
2726 xmatri_offset = xmatri_dim1 + 1;
2727 xmatri -= xmatri_offset;
2728 tabaux_dim1 = *nderiv + 1 + *ndimen;
2729 tabaux_offset = tabaux_dim1 + 1;
2730 tabaux -= tabaux_offset;
2731 ctrtes_dim1 = *ndimen;
2732 ctrtes_offset = ctrtes_dim1 * 3 + 1;
2733 ctrtes -= ctrtes_offset;
2736 ibb = AdvApp2Var_SysBase::mnfndeb_();
2738 AdvApp2Var_SysBase::mgenmsg_("MMCVCTX", 7L);
2741 AdvApp2Var_MathBase::mmeps1_(&eps1);
2743 /* ****************** CALCULATION OF EVEN COEFFICIENTS *********************
2745 /* ------------------------- Initialization -----------------------------
2748 nordr = *nderiv + 1;
2750 for (ncf = 1; ncf <= i__1; ++ncf) {
2751 tabaux[ncf + tabaux_dim1] = 1.;
2755 /* ---------------- Calculation of terms corresponding to derivatives -------
2759 for (ndv = 2; ndv <= i__1; ++ndv) {
2761 for (ncf = 1; ncf <= i__2; ++ncf) {
2762 tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
2763 tabaux_dim1] * (doublereal) ((ncf << 1) - ndv);
2769 /* ------------------ Writing the second member -----------------------
2774 for (ndv = 1; ndv <= i__1; ++ndv) {
2776 for (nd = 1; nd <= i__2; ++nd) {
2777 tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
2778 + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2779 * ctrtes_dim1]) / 2.;
2786 /* -------------------- Resolution of the system ---------------------------
2789 mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2790 xmatri_offset], iercod);
2795 for (nd = 1; nd <= i__1; ++nd) {
2797 for (ncf = 1; ncf <= i__2; ++ncf) {
2798 crvres[(ncf << 1) - 1 + nd * crvres_dim1] = xmatri[ncf + nd *
2805 /* ***************** CALCULATION OF UNEVEN COEFFICIENTS ********************
2807 /* ------------------------- Initialization -----------------------------
2812 for (ncf = 1; ncf <= i__1; ++ncf) {
2813 tabaux[ncf + tabaux_dim1] = 1.;
2817 /* ---------------- Calculation of terms corresponding to derivatives -------
2821 for (ndv = 2; ndv <= i__1; ++ndv) {
2823 for (ncf = 1; ncf <= i__2; ++ncf) {
2824 tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
2825 tabaux_dim1] * (doublereal) ((ncf << 1) - ndv + 1);
2831 /* ------------------ Writing of the second member -----------------------
2836 for (ndv = 1; ndv <= i__1; ++ndv) {
2838 for (nd = 1; nd <= i__2; ++nd) {
2839 tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
2840 + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2841 * ctrtes_dim1]) / 2.;
2848 /* -------------------- Solution of the system ---------------------------
2851 mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2852 xmatri_offset], iercod);
2857 for (nd = 1; nd <= i__1; ++nd) {
2859 for (ncf = 1; ncf <= i__2; ++ncf) {
2860 crvres[(ncf << 1) + nd * crvres_dim1] = xmatri[ncf + nd *
2867 /* --------------------------- The end ----------------------------------
2872 AdvApp2Var_SysBase::maermsg_("MMCVCTX", iercod, 7L);
2875 AdvApp2Var_SysBase::mgsomsg_("MMCVCTX", 7L);
2881 //=======================================================================
2882 //function : AdvApp2Var_MathBase::mmcvinv_
2884 //=======================================================================
2885 int AdvApp2Var_MathBase::mmcvinv_(integer *ndimax,
2892 /* Initialized data */
2894 static char nomprg[8+1] = "MMCVINV ";
2896 /* System generated locals */
2897 integer curve_dim1, curve_offset, curveo_dim1, curveo_offset, i__1, i__2;
2899 /* Local variables */
2900 integer i__, nd, ibb;
2903 /* ***********************************************************************
2908 /* Inversion of arguments of the final curve. */
2912 /* SMOOTHING,CURVE */
2915 /* INPUT ARGUMENTS : */
2916 /* ------------------ */
2918 /* NDIM: Space Dimension. */
2919 /* NCOEF: Degree of the polynom. */
2920 /* CURVEO: The curve before inversion. */
2922 /* OUTPUT ARGUMENTS : */
2923 /* ------------------- */
2924 /* CURVE: The curve after inversion. */
2926 /* COMMONS USED : */
2927 /* ---------------- */
2928 /* REFERENCES APPELEES : */
2929 /* ----------------------- */
2930 /* DESCRIPTION/NOTES/LIMITATIONS : */
2931 /* ----------------------------------- */
2932 /* ***********************************************************************
2935 /* The name of the routine */
2936 /* Parameter adjustments */
2937 curve_dim1 = *ndimax;
2938 curve_offset = curve_dim1 + 1;
2939 curve -= curve_offset;
2940 curveo_dim1 = *ncoef;
2941 curveo_offset = curveo_dim1 + 1;
2942 curveo -= curveo_offset;
2946 ibb = AdvApp2Var_SysBase::mnfndeb_();
2948 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
2952 for (i__ = 1; i__ <= i__1; ++i__) {
2954 for (nd = 1; nd <= i__2; ++nd) {
2955 curve[nd + i__ * curve_dim1] = curveo[i__ + nd * curveo_dim1];
2964 //=======================================================================
2965 //function : mmcvstd_
2967 //=======================================================================
2968 int mmcvstd_(integer *ncofmx,
2976 /* System generated locals */
2977 integer courbe_dim1, crvcan_dim1, crvcan_offset, i__1, i__2, i__3;
2979 /* Local variables */
2980 integer ndeg, i__, j, j1, nd, ibb;
2984 /* ***********************************************************************
2989 /* Transform curve defined between [-1,1] into [0,1]. */
2993 /* LIMITATION,RESTRICTION,CURVE */
2995 /* INPUT ARGUMENTS : */
2996 /* ------------------ */
2997 /* NDIMAX : Dimension of the space. */
2998 /* NDIMEN : Dimension of the curve. */
2999 /* NCOEFF : Degree of the curve. */
3000 /* CRVCAN(NCOFMX,NDIMEN): The curve is defined at the interval [-1,1]. */
3002 /* OUTPUT ARGUMENTS : */
3003 /* ------------------- */
3004 /* CURVE(NDIMAX,NCOEFF): Curve defined at the interval [0,1]. */
3006 /* COMMONS USED : */
3007 /* ---------------- */
3009 /* REFERENCES CALLED : */
3010 /* ----------------------- */
3012 /* DESCRIPTION/NOTES/LIMITATIONS : */
3013 /* ----------------------------------- */
3015 /* ***********************************************************************
3018 /* Name of the program. */
3021 /* **********************************************************************
3026 /* Provides binomial coefficients (Pascal triangle). */
3030 /* Binomial coefficient from 0 to 60. read only . init by block data */
3032 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3033 /* ----------------------------------- */
3034 /* Binomial coefficients form a triangular matrix. */
3035 /* This matrix is completed in table CNP by its transposition. */
3036 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
3038 /* Initialization is done with block-data MMLLL09.RES, */
3039 /* created by the program MQINICNP.FOR. */
3041 /* **********************************************************************
3046 /* ***********************************************************************
3049 /* Parameter adjustments */
3050 courbe_dim1 = *ndimax;
3052 crvcan_dim1 = *ncofmx;
3053 crvcan_offset = crvcan_dim1;
3054 crvcan -= crvcan_offset;
3057 ibb = AdvApp2Var_SysBase::mnfndeb_();
3059 AdvApp2Var_SysBase::mgenmsg_("MMCVSTD", 7L);
3063 /* ------------------ Construction of the resulting curve ----------------
3067 for (nd = 1; nd <= i__1; ++nd) {
3069 for (j = 0; j <= i__2; ++j) {
3072 for (i__ = j; i__ <= i__3; i__ += 2) {
3073 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
3077 courbe[nd + j * courbe_dim1] = bid;
3082 for (i__ = j1; i__ <= i__3; i__ += 2) {
3083 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
3087 courbe[nd + j * courbe_dim1] -= bid;
3093 /* ------------------- Renormalization of the CURVE -------------------------
3098 for (i__ = 0; i__ <= i__1; ++i__) {
3100 for (nd = 1; nd <= i__2; ++nd) {
3101 courbe[nd + i__ * courbe_dim1] *= bid;
3108 /* ----------------------------- The end --------------------------------
3112 AdvApp2Var_SysBase::mgsomsg_("MMCVSTD", 7L);
3117 //=======================================================================
3118 //function : AdvApp2Var_MathBase::mmdrc11_
3120 //=======================================================================
3121 int AdvApp2Var_MathBase::mmdrc11_(integer *iordre,
3126 doublereal *mfactab)
3129 /* System generated locals */
3130 integer courbe_dim1, courbe_offset, points_dim2, points_offset, i__1,
3133 /* Local variables */
3135 integer ndeg, i__, j, ndgcb, nd, ibb;
3138 /* **********************************************************************
3143 /* Calculation of successive derivatives of equation CURVE with */
3144 /* parameters -1, 1 from order 0 to order IORDRE */
3145 /* included. The calculation is produced without knowing the coefficients of */
3146 /* derivatives of the curve. */
3150 /* POSITIONING,EXTREMITIES,CURVE,DERIVATIVE. */
3152 /* INPUT ARGUMENTS : */
3153 /* ------------------ */
3154 /* IORDRE : Maximum order of calculation of derivatives. */
3155 /* NDIMEN : Dimension of the space. */
3156 /* NCOEFF : Number of coefficients of the curve (degree+1). */
3157 /* COURBE : Table of coefficients of the curve. */
3159 /* OUTPUT ARGUMENTS : */
3160 /* ------------------- */
3161 /* POINTS : Table of values of consecutive derivatives */
3162 /* of parameters -1.D0 and 1.D0. */
3163 /* MFACTAB : Auxiliary table for calculation of factorial(I).
3166 /* COMMONS USED : */
3167 /* ---------------- */
3170 /* REFERENCES CALLED : */
3171 /* ----------------------- */
3173 /* DESCRIPTION/NOTES/LIMITATIONS : */
3174 /* ----------------------------------- */
3176 /* ---> ATTENTION, the coefficients of the curve are */
3177 /* in a reverse order. */
3179 /* ---> The algorithm of calculation of derivatives is based on */
3180 /* generalization of Horner scheme : */
3182 /* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
3185 /* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
3187 /* aj = a(j-1).x + u(k-j) */
3188 /* bj = b(j-1).x + a(j-1) */
3189 /* cj = c(j-1).x + b(j-1) */
3191 /* So : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
3193 /* The algorithm is generalized easily for calculation of */
3200 /* Reference : D. KNUTH, "The Art of Computer Programming" */
3201 /* --------- Vol. 2/Seminumerical Algorithms */
3202 /* Addison-Wesley Pub. Co. (1969) */
3203 /* pages 423-425. */
3205 /* **********************************************************************
3208 /* Name of the routine */
3210 /* Parameter adjustments */
3211 points_dim2 = *iordre + 1;
3212 points_offset = (points_dim2 << 1) + 1;
3213 points -= points_offset;
3214 courbe_dim1 = *ncoeff;
3215 courbe_offset = courbe_dim1;
3216 courbe -= courbe_offset;
3219 ibb = AdvApp2Var_SysBase::mnfndeb_();
3221 AdvApp2Var_SysBase::mgenmsg_("MMDRC11", 7L);
3224 if (*iordre < 0 || *ncoeff < 1) {
3228 /* ------------------- Initialization of table POINTS -----------------
3231 ndgcb = *ncoeff - 1;
3233 for (nd = 1; nd <= i__1; ++nd) {
3234 points[(nd * points_dim2 << 1) + 1] = courbe[ndgcb + nd * courbe_dim1]
3236 points[(nd * points_dim2 << 1) + 2] = courbe[ndgcb + nd * courbe_dim1]
3242 for (nd = 1; nd <= i__1; ++nd) {
3244 for (j = 1; j <= i__2; ++j) {
3245 points[((j + nd * points_dim2) << 1) + 1] = 0.;
3246 points[((j + nd * points_dim2) << 1) + 2] = 0.;
3252 /* Calculation with parameter -1 and 1 */
3255 for (nd = 1; nd <= i__1; ++nd) {
3257 for (ndeg = 1; ndeg <= i__2; ++ndeg) {
3258 for (i__ = *iordre; i__ >= 1; --i__) {
3259 points[((i__ + nd * points_dim2) << 1) + 1] = -points[((i__ + nd
3260 * points_dim2) << 1) + 1] + points[((i__ - 1 + nd *
3261 points_dim2) << 1) + 1];
3262 points[((i__ + nd * points_dim2) << 1) + 2] += points[((i__ - 1
3263 + nd * points_dim2) << 1) + 2];
3266 points[(nd * points_dim2 << 1) + 1] = -points[(nd * points_dim2 <<
3267 1) + 1] + courbe[ndgcb - ndeg + nd * courbe_dim1];
3268 points[(nd * points_dim2 << 1) + 2] += courbe[ndgcb - ndeg + nd *
3275 /* --------------------- Multiplication by factorial(I) --------------
3279 mfac_(&mfactab[1], iordre);
3282 for (nd = 1; nd <= i__1; ++nd) {
3284 for (i__ = 2; i__ <= i__2; ++i__) {
3285 points[((i__ + nd * points_dim2) << 1) + 1] = mfactab[i__] *
3286 points[((i__ + nd * points_dim2) << 1) + 1];
3287 points[((i__ + nd * points_dim2) << 1) + 2] = mfactab[i__] *
3288 points[((i__ + nd * points_dim2) << 1) + 2];
3295 /* ---------------------------- End -------------------------------------
3300 AdvApp2Var_SysBase::mgsomsg_("MMDRC11", 7L);
3305 //=======================================================================
3306 //function : mmdrvcb_
3308 //=======================================================================
3309 int mmdrvcb_(integer *ideriv,
3318 /* System generated locals */
3319 integer courbe_dim1, tabpnt_dim1, i__1, i__2, i__3;
3321 /* Local variables */
3322 integer ndeg, i__, j, nd, ndgcrb, iptpnt, ibb;
3325 /* ***********************************************************************
3329 /* Calculation of successive derivatives of equation CURVE with */
3330 /* parameter TPARAM from order 0 to order IDERIV included. */
3331 /* The calculation is produced without knowing the coefficients of */
3332 /* derivatives of the CURVE. */
3336 /* POSITIONING,PARAMETER,CURVE,DERIVATIVE. */
3338 /* INPUT ARGUMENTS : */
3339 /* ------------------ */
3340 /* IORDRE : Maximum order of calculation of derivatives. */
3341 /* NDIMEN : Dimension of the space. */
3342 /* NCOEFF : Number of coefficients of the curve (degree+1). */
3343 /* COURBE : Table of coefficients of the curve. */
3344 /* TPARAM : Value of the parameter where the curve should be evaluated. */
3346 /* OUTPUT ARGUMENTS : */
3347 /* ------------------- */
3348 /* TABPNT : Table of values of consecutive derivatives */
3349 /* of parameter TPARAM. */
3350 /* IERCOD : 0 = OK, */
3351 /* 1 = incoherent input. */
3353 /* COMMONS USED : */
3354 /* ---------------- */
3357 /* REFERENCES CALLED : */
3358 /* ----------------------- */
3360 /* DESCRIPTION/NOTES/LIMITATIONS : */
3361 /* ----------------------------------- */
3363 /* The algorithm of calculation of derivatives is based on */
3364 /* generalization of the Horner scheme : */
3366 /* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
3369 /* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
3371 /* aj = a(j-1).x + u(k-j) */
3372 /* bj = b(j-1).x + a(j-1) */
3373 /* cj = c(j-1).x + b(j-1) */
3375 /* So, it is obtained : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
3377 /* The algorithm can be easily generalized for the calculation of */
3384 /* Reference : D. KNUTH, "The Art of Computer Programming" */
3385 /* --------- Vol. 2/Seminumerical Algorithms */
3386 /* Addison-Wesley Pub. Co. (1969) */
3387 /* pages 423-425. */
3389 /* ---> To evaluare derivatives at 0 and 1, it is preferable */
3390 /* to use routine MDRV01.FOR . */
3392 /* **********************************************************************
3395 /* Name of the routine */
3397 /* Parameter adjustments */
3398 tabpnt_dim1 = *ndim;
3400 courbe_dim1 = *ndim;
3404 ibb = AdvApp2Var_SysBase::mnfndeb_();
3406 AdvApp2Var_SysBase::mgenmsg_("MMDRVCB", 7L);
3409 if (*ideriv < 0 || *ncoeff < 1) {
3415 /* ------------------- Initialization of table TABPNT -----------------
3418 ndgcrb = *ncoeff - 1;
3420 for (nd = 1; nd <= i__1; ++nd) {
3421 tabpnt[nd] = courbe[nd + ndgcrb * courbe_dim1];
3428 iptpnt = *ndim * *ideriv;
3429 AdvApp2Var_SysBase::mvriraz_(&iptpnt,
3430 &tabpnt[tabpnt_dim1 + 1]);
3433 /* ------------------------ Calculation of parameter TPARAM ------------------
3437 for (ndeg = 1; ndeg <= i__1; ++ndeg) {
3439 for (nd = 1; nd <= i__2; ++nd) {
3440 for (i__ = *ideriv; i__ >= 1; --i__) {
3441 tabpnt[nd + i__ * tabpnt_dim1] = tabpnt[nd + i__ *
3442 tabpnt_dim1] * *tparam + tabpnt[nd + (i__ - 1) *
3446 tabpnt[nd] = tabpnt[nd] * *tparam + courbe[nd + (ndgcrb - ndeg) *
3453 /* --------------------- Multiplication by factorial(I) -------------
3457 for (i__ = 2; i__ <= i__1; ++i__) {
3459 for (j = 2; j <= i__2; ++j) {
3461 for (nd = 1; nd <= i__3; ++nd) {
3462 tabpnt[nd + i__ * tabpnt_dim1] = (doublereal) j * tabpnt[nd +
3471 /* --------------------------- The end ---------------------------------
3476 AdvApp2Var_SysBase::maermsg_("MMDRVCB", iercod, 7L);
3481 //=======================================================================
3482 //function : AdvApp2Var_MathBase::mmdrvck_
3484 //=======================================================================
3485 int AdvApp2Var_MathBase::mmdrvck_(integer *ncoeff,
3493 /* Initialized data */
3495 static doublereal mmfack[21] = { 1.,2.,6.,24.,120.,720.,5040.,40320.,
3496 362880.,3628800.,39916800.,479001600.,6227020800.,87178291200.,
3497 1.307674368e12,2.0922789888e13,3.55687428096e14,6.402373705728e15,
3498 1.21645100408832e17,2.43290200817664e18,5.109094217170944e19 };
3500 /* System generated locals */
3501 integer courbe_dim1, courbe_offset, i__1, i__2;
3503 /* Local variables */
3504 integer i__, j, k, nd;
3505 doublereal mfactk, bid;
3508 /* IMPLICIT INTEGER (I-N) */
3509 /* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
3512 /* ***********************************************************************
3517 /* Calculate the value of a derived curve of order IDERIV in */
3518 /* a point of parameter TPARAM. */
3522 /* POSITIONING,CURVE,DERIVATIVE of ORDER K. */
3524 /* INPUT ARGUMENTS : */
3525 /* ------------------ */
3526 /* NCOEFF : Degree +1 of the curve. */
3527 /* NDIMEN : Dimension of the space (2 or 3 in general) */
3528 /* COURBE : Table of coefficients of the curve. */
3529 /* IDERIV : Required order of derivation : 1=1st derivative, etc... */
3530 /* TPARAM : Value of parameter of the curve. */
3532 /* OUTPUT ARGUMENTS : */
3533 /* ------------------- */
3534 /* PNTCRB : Point of parameter TPARAM on the derivative of order */
3535 /* IDERIV of CURVE. */
3537 /* COMMONS USED : */
3538 /* ---------------- */
3541 /* REFERENCES CALLED : */
3542 /* ---------------------- */
3544 /* DESCRIPTION/NOTES/LIMITATIONS : */
3545 /* ----------------------------------- */
3547 /* The code below was written basing on the following algorithm :
3550 /* Let P(t) = a1 + a2*t + ... an*t**n. The derivative of order k of P */
3551 /* (containing n-k coefficients) is calculated as follows : */
3553 /* Pk(t) = a(k+1)*CNP(k,k)*k! */
3554 /* + a(k+2)*CNP(k+1,k)*k! * t */
3558 /* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
3560 /* Evaluation is produced following the classic Horner scheme. */
3562 /* ***********************************************************************
3566 /* Factorials (1 to 21) caculated on VAX in R*16 */
3569 /* **********************************************************************
3574 /* Serves to provide binomial coefficients (Pascal triangle). */
3578 /* Binomial Coeff from 0 to 60. read only . init by block data */
3580 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3581 /* ----------------------------------- */
3582 /* Binomial coefficients form a triangular matrix. */
3583 /* This matrix is completed in table CNP by its transposition. */
3584 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
3586 /* Initialization is done by block-data MMLLL09.RES, */
3587 /* created by program MQINICNP.FOR. */
3589 /* **********************************************************************
3594 /* ***********************************************************************
3597 /* Parameter adjustments */
3599 courbe_dim1 = *ndimen;
3600 courbe_offset = courbe_dim1 + 1;
3601 courbe -= courbe_offset;
3605 /* -------------- Case when the order of derivative is greater than -------------------
3607 /* ---------------- the degree of the curve ---------------------
3610 if (*ideriv >= *ncoeff) {
3612 for (nd = 1; nd <= i__1; ++nd) {
3618 /* **********************************************************************
3620 /* General processing*/
3621 /* **********************************************************************
3623 /* --------------------- Calculation of Factorial(IDERIV) ------------------
3627 if (*ideriv <= 21 && *ideriv > 0) {
3628 mfactk = mmfack[k - 1];
3632 for (i__ = 2; i__ <= i__1; ++i__) {
3638 /* ------- Calculation of derivative of order IDERIV of CURVE in TPARAM -----
3640 /* ---> Attention : binomial coefficient C(n,m) is represented in */
3641 /* MCCNP by CNP(N,M). */
3644 for (nd = 1; nd <= i__1; ++nd) {
3645 pntcrb[nd] = courbe[nd + *ncoeff * courbe_dim1] * mmcmcnp_.cnp[*
3646 ncoeff - 1 + k * 61] * mfactk;
3651 for (j = *ncoeff - 1; j >= i__1; --j) {
3652 bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
3654 for (nd = 1; nd <= i__2; ++nd) {
3655 pntcrb[nd] = pntcrb[nd] * *tparam + courbe[nd + j * courbe_dim1] *
3662 /* -------------------------------- The end -----------------------------
3670 //=======================================================================
3671 //function : AdvApp2Var_MathBase::mmeps1_
3673 //=======================================================================
3674 int AdvApp2Var_MathBase::mmeps1_(doublereal *epsilo)
3677 /* ***********************************************************************
3682 /* Extraction of EPS1 from COMMON MPRCSN. EPS1 is spatial zero */
3683 /* equal to 1.D-9 */
3687 /* MPRCSN,PRECISON,EPS1. */
3689 /* INPUT ARGUMENTS : */
3690 /* ------------------ */
3693 /* OUTPUT ARGUMENTS : */
3694 /* ------------------- */
3695 /* EPSILO : Value of EPS1 (spatial zero (10**-9)) */
3697 /* COMMONS USED : */
3698 /* ---------------- */
3700 /* REFERENCES CALLED : */
3701 /* ----------------------- */
3703 /* DESCRIPTION/NOTES/LIMITATIONS : */
3704 /* ----------------------------------- */
3705 /* EPS1 is ABSOLUTE spatial zero, so it is necessary */
3706 /* to use it whenever it is necessary to test if a variable */
3707 /* is null. For example, if the norm of a vector is lower than */
3708 /* EPS1, this vector is NULL ! (when one works in */
3709 /* REAL*8) It is absolutely not advised to test arguments */
3710 /* compared to EPS1**2. Taking into account the rounding errors inevitable */
3711 /* during calculations, this causes testing compared to 0.D0. */
3713 /* ***********************************************************************
3718 /* ***********************************************************************
3723 /* Gives tolerances of invalidity in stream */
3724 /* as well as limits of iterative processes */
3726 /* general context, modifiable by the user */
3730 /* PARAMETER , TOLERANCE */
3732 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3733 /* ----------------------------------- */
3734 /* INITIALISATION : profile , **VIA MPRFTX** at input in stream
3735 /* loading of default values of the profile in MPRFTX at input */
3736 /* in stream. They are preserved in local variables of MPRFTX */
3738 /* Reset of default values : MDFINT */
3739 /* Interactive modification by the user : MDBINT */
3741 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
3742 /* MEPSPB ... EPS3,EPS4 */
3743 /* MEPSLN ... EPS2, NITERM , NITERR */
3744 /* MEPSNR ... EPS2 , NITERM */
3745 /* MITERR ... NITERR */
3747 /* ***********************************************************************
3750 /* NITERM : max nb of iterations */
3751 /* NITERR : nb of rapid iterations */
3752 /* EPS1 : tolerance of 3D null distance */
3753 /* EPS2 : tolerance of parametric null distance */
3754 /* EPS3 : tolerance to avoid division by 0.. */
3755 /* EPS4 : angular tolerance */
3759 /* ***********************************************************************
3761 *epsilo = mmprcsn_.eps1;
3766 //=======================================================================
3767 //function : mmexthi_
3769 //=======================================================================
3770 int mmexthi_(integer *ndegre,
3774 /* System generated locals */
3777 /* Local variables */
3778 integer iadd, ideb, ndeg2, nmod2, ii, ibb;
3781 /* **********************************************************************
3786 /* Extract of common LDGRTL the weight of formulas of */
3787 /* Gauss quadrature on all roots of Legendre polynoms of degree */
3788 /* NDEGRE defined on [-1,1]. */
3792 /* ALL, AB_SPECIFI::COMMON&, EXTRACTION, &WEIGHT, &GAUSS. */
3794 /* INPUT ARGUMENTS : */
3795 /* ------------------ */
3796 /* NDEGRE : Mathematic degree of Legendre polynom. It should have */
3797 /* 2 <= NDEGRE <= 61. */
3799 /* OUTPUT ARGUMENTS : */
3800 /* ------------------- */
3801 /* HWGAUS : The table of weights of Gauss quadrature formulas */
3802 /* relative to NDEGRE roots of a polynome de Legendre de */
3805 /* COMMONS UTILISES : */
3806 /* ---------------- */
3809 /* REFERENCES CALLED : */
3810 /* ----------------------- */
3812 /* DESCRIPTION/NOTES/LIMITATIONS : */
3813 /* ----------------------------------- */
3814 /* ATTENTION: The condition on NDEGRE ( 2 <= NDEGRE <= 61) is not */
3815 /* tested. The caller should make the test.
3817 /* Name of the routine */
3820 /* Common MLGDRTL: */
3821 /* This common includes POSITIVE roots of Legendre polynims */
3822 /* AND weights of Gauss quadrature formulas on all */
3823 /* POSITIVE roots of Legendre polynoms. */
3827 /* ***********************************************************************
3832 /* The common of Legendre roots. */
3838 /* DESCRIPTION/NOTES/LIMITATIONS : */
3839 /* ----------------------------------- */
3841 /* ***********************************************************************
3847 /* ROOTAB : Table of all roots of Legendre polynoms */
3848 /* within the interval [0,1]. They are ranked for the degrees increasing from */
3850 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
3851 /* The adressing is the same. */
3852 /* HI0TAB : Table of Legendre interpolators for root x=0 */
3853 /* of polynoms of UNEVEN degree. */
3854 /* RTLTB0 : Table of Li(uk) where uk are the roots of */
3855 /* Legendre polynom of EVEN degree. */
3856 /* RTLTB1 : Table of Li(uk) where uk are the roots of */
3857 /* Legendre polynom of UNEVEN degree. */
3860 /************************************************************************
3862 /* Parameter adjustments */
3866 ibb = AdvApp2Var_SysBase::mnfndeb_();
3868 AdvApp2Var_SysBase::mgenmsg_("MMEXTHI", 7L);
3871 ndeg2 = *ndegre / 2;
3872 nmod2 = *ndegre % 2;
3874 /* Address of Gauss weight associated to the 1st strictly */
3875 /* positive root of Legendre polynom of degree NDEGRE in MLGDRTL. */
3877 iadd = ndeg2 * (ndeg2 - 1) / 2 + 1;
3879 /* Index of the 1st HWGAUS element associated to the 1st strictly */
3880 /* positive root of Legendre polynom of degree NDEGRE. */
3882 ideb = (*ndegre + 1) / 2 + 1;
3884 /* Reading of weights associated to strictly positive roots. */
3887 for (ii = ideb; ii <= i__1; ++ii) {
3888 kpt = iadd + ii - ideb;
3889 hwgaus[ii] = mlgdrtl_.hiltab[kpt + nmod2 * 465 - 1];
3893 /* For strictly negative roots, the weight is the same. */
3894 /* i.e HW(1) = HW(NDEGRE), HW(2) = HW(NDEGRE-1), etc... */
3897 for (ii = 1; ii <= i__1; ++ii) {
3898 hwgaus[ii] = hwgaus[*ndegre + 1 - ii];
3902 /* Case of uneven NDEGRE, 0 is root of Legendre polynom, */
3903 /* associated Gauss weights are loaded. */
3906 hwgaus[ndeg2 + 1] = mlgdrtl_.hi0tab[ndeg2];
3909 /* --------------------------- The end ----------------------------------
3913 AdvApp2Var_SysBase::mgsomsg_("MMEXTHI", 7L);
3918 //=======================================================================
3919 //function : mmextrl_
3921 //=======================================================================
3922 int mmextrl_(integer *ndegre,
3925 /* System generated locals */
3928 /* Local variables */
3929 integer iadd, ideb, ndeg2, nmod2, ii, ibb;
3933 /* **********************************************************************
3938 /* Extract of the Common LDGRTL of Legendre polynom roots */
3939 /* of degree NDEGRE defined on [-1,1]. */
3943 /* ALL, AB_SPECIFI::COMMON&, EXTRACTION, &ROOT, &LEGENDRE. */
3945 /* INPUT ARGUMENTS : */
3946 /* ------------------ */
3947 /* NDEGRE : Mathematic degree of Legendre polynom. */
3948 /* It is required to have 2 <= NDEGRE <= 61. */
3950 /* OUTPUT ARGUMENTS : */
3951 /* ------------------- */
3952 /* ROOTLG : The table of roots of Legendre polynom of degree */
3953 /* NDEGRE defined on [-1,1]. */
3955 /* COMMONS USED : */
3956 /* ---------------- */
3959 /* REFERENCES CALLED : */
3960 /* ----------------------- */
3962 /* DESCRIPTION/NOTES/LIMITATIONS : */
3963 /* ----------------------------------- */
3964 /* ATTENTION: Condition of NDEGRE ( 2 <= NDEGRE <= 61) is not */
3965 /* tested. The caller should make the test. */
3967 /* **********************************************************************
3971 /* Name of the routine */
3974 /* Common MLGDRTL: */
3975 /* This common includes POSITIVE roots of Legendre polynoms */
3976 /* AND the weight of Gauss quadrature formulas on all */
3977 /* POSITIVE roots of Legendre polynoms. */
3979 /* ***********************************************************************
3984 /* The common of Legendre roots. */
3991 /* ***********************************************************************
3994 /* ROOTAB : Table of all roots of Legendre polynoms */
3995 /* within the interval [0,1]. They are ranked for the degrees increasing from */
3997 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
3998 /* The adressing is the same. */
3999 /* HI0TAB : Table of Legendre interpolators for root x=0 */
4000 /* of polynoms of UNEVEN degree. */
4001 /* RTLTB0 : Table of Li(uk) where uk are the roots of */
4002 /* Legendre polynom of EVEN degree. */
4003 /* RTLTB1 : Table of Li(uk) where uk are the roots of */
4004 /* Legendre polynom of UNEVEN degree. */
4007 /************************************************************************
4009 /* Parameter adjustments */
4013 ibb = AdvApp2Var_SysBase::mnfndeb_();
4015 AdvApp2Var_SysBase::mgenmsg_("MMEXTRL", 7L);
4018 ndeg2 = *ndegre / 2;
4019 nmod2 = *ndegre % 2;
4021 /* Address of the 1st strictly positive root of Legendre polynom */
4022 /* of degree NDEGRE in MLGDRTL. */
4024 iadd = ndeg2 * (ndeg2 - 1) / 2 + 1;
4026 /* Indice, in ROOTLG, of the 1st strictly positive root */
4027 /* of Legendre polynom of degree NDEGRE. */
4029 ideb = (*ndegre + 1) / 2 + 1;
4031 /* Reading of strictly positive roots. */
4034 for (ii = ideb; ii <= i__1; ++ii) {
4035 kpt = iadd + ii - ideb;
4036 rootlg[ii] = mlgdrtl_.rootab[kpt + nmod2 * 465 - 1];
4040 /* Strictly negative roots are equal to positive roots
4042 /* to the sign i.e RT(1) = -RT(NDEGRE), RT(2) = -RT(NDEGRE-1), etc...
4046 for (ii = 1; ii <= i__1; ++ii) {
4047 rootlg[ii] = -rootlg[*ndegre + 1 - ii];
4051 /* Case NDEGRE uneven, 0 is root of Legendre polynom. */
4054 rootlg[ndeg2 + 1] = 0.;
4057 /* -------------------------------- THE END -----------------------------
4061 AdvApp2Var_SysBase::mgenmsg_("MMEXTRL", 7L);
4066 //=======================================================================
4067 //function : AdvApp2Var_MathBase::mmfmca8_
4069 //=======================================================================
4070 int AdvApp2Var_MathBase::mmfmca8_(const integer *ndimen,
4071 const integer *ncoefu,
4072 const integer *ncoefv,
4073 const integer *ndimax,
4074 const integer *ncfumx,
4075 const integer *,//ncfvmx,
4080 /* System generated locals */
4081 integer tabini_dim1, tabini_dim2, tabini_offset, tabres_dim1, tabres_dim2,
4084 /* Local variables */
4085 integer i__, j, k, ilong;
4089 /* **********************************************************************
4094 /* Expansion of a table containing only most important things into a */
4095 /* greater data table. */
4099 /* ALL, MATH_ACCES:: CARREAU&, DECOMPRESSION, &CARREAU */
4101 /* INPUT ARGUMENTS : */
4102 /* ------------------ */
4103 /* NDIMEN: Dimension of the workspace. */
4104 /* NCOEFU: Degree +1 of the table by u. */
4105 /* NCOEFV: Degree +1 of the table by v. */
4106 /* NDIMAX: Max dimension of the space. */
4107 /* NCFUMX: Max Degree +1 of the table by u. */
4108 /* NCFVMX: Max Degree +1 of the table by v. */
4109 /* TABINI: The table to be decompressed. */
4111 /* OUTPUT ARGUMENTS : */
4112 /* ------------------- */
4113 /* TABRES: Decompressed table. */
4115 /* COMMONS USED : */
4116 /* ---------------- */
4118 /* REFERENCES CALLED : */
4119 /* ----------------------- */
4121 /* DESCRIPTION/NOTES/LIMITATIONS : */
4122 /* ----------------------------------- */
4123 /* The following call : */
4125 /* CALL MMFMCA8(NDIMEN,NCOEFU,NCOEFV,NDIMAX,NCFUMX,NCFVMX,TABINI,TABINI)
4128 /* where TABINI is input/output argument, is possible provided */
4129 /* that the caller has declared TABINI in (NDIMAX,NCFUMX,NCFVMX) */
4131 /* ATTENTION : it is not checked that NDIMAX >= NDIMEN, */
4132 /* NCOEFU >= NCFMXU and NCOEFV >= NCFMXV. */
4134 /* **********************************************************************
4138 /* Parameter adjustments */
4139 tabini_dim1 = *ndimen;
4140 tabini_dim2 = *ncoefu;
4141 tabini_offset = tabini_dim1 * (tabini_dim2 + 1) + 1;
4142 tabini -= tabini_offset;
4143 tabres_dim1 = *ndimax;
4144 tabres_dim2 = *ncfumx;
4145 tabres_offset = tabres_dim1 * (tabres_dim2 + 1) + 1;
4146 tabres -= tabres_offset;
4149 if (*ndimax == *ndimen) {
4153 /* ----------------------- decompression NDIMAX<>NDIMEN -----------------
4156 for (k = *ncoefv; k >= 1; --k) {
4157 for (j = *ncoefu; j >= 1; --j) {
4158 for (i__ = *ndimen; i__ >= 1; --i__) {
4159 tabres[i__ + (j + k * tabres_dim2) * tabres_dim1] = tabini[
4160 i__ + (j + k * tabini_dim2) * tabini_dim1];
4169 /* ----------------------- decompression NDIMAX=NDIMEN ------------------
4173 if (*ncoefu == *ncfumx) {
4176 ilong = (*ndimen << 3) * *ncoefu;
4177 for (k = *ncoefv; k >= 1; --k) {
4178 AdvApp2Var_SysBase::mcrfill_(&ilong,
4179 &tabini[(k * tabini_dim2 + 1) * tabini_dim1 + 1],
4180 &tabres[(k * tabres_dim2 + 1) * tabres_dim1 + 1]);
4185 /* ----------------- decompression NDIMAX=NDIMEN,NCOEFU=NCFUMX ----------
4189 ilong = (*ndimen << 3) * *ncoefu * *ncoefv;
4190 AdvApp2Var_SysBase::mcrfill_(&ilong,
4191 &tabini[tabini_offset],
4192 &tabres[tabres_offset]);
4195 /* ---------------------------- The end ---------------------------------
4202 //=======================================================================
4203 //function : AdvApp2Var_MathBase::mmfmca9_
4205 //=======================================================================
4206 int AdvApp2Var_MathBase::mmfmca9_(integer *ndimax,
4216 /* System generated locals */
4217 integer tabini_dim1, tabini_dim2, tabini_offset, tabres_dim1, tabres_dim2,
4218 tabres_offset, i__1, i__2, i__3;
4220 /* Local variables */
4221 integer i__, j, k, ilong;
4225 /* **********************************************************************
4230 /* Compression of a data table in a table */
4231 /* containing only the main data (the input table is not removed). */
4235 /* ALL, MATH_ACCES:: CARREAU&, COMPRESSION, &CARREAU */
4237 /* INPUT ARGUMENTS : */
4238 /* ------------------ */
4239 /* NDIMAX: Max dimension of the space. */
4240 /* NCFUMX: Max degree +1 of the table by u. */
4241 /* NCFVMX: Max degree +1 of the table by v. */
4242 /* NDIMEN: Dimension of the workspace. */
4243 /* NCOEFU: Degree +1 of the table by u. */
4244 /* NCOEFV: Degree +1 of the table by v. */
4245 /* TABINI: The table to compress. */
4247 /* OUTPUT ARGUMENTS : */
4248 /* ------------------- */
4249 /* TABRES: The compressed table. */
4251 /* COMMONS USED : */
4252 /* ---------------- */
4254 /* REFERENCES CALLED : */
4255 /* ----------------------- */
4257 /* DESCRIPTION/NOTES/LIMITATIONS : */
4258 /* ----------------------------------- */
4259 /* The following call : */
4261 /* CALL MMFMCA9(NDIMAX,NCFUMX,NCFVMX,NDIMEN,NCOEFU,NCOEFV,TABINI,TABINI)
4264 /* where TABINI is input/output argument, is possible provided */
4265 /* that the caller has checked that : */
4267 /* NDIMAX > NDIMEN, */
4268 /* or NDIMAX = NDIMEN and NCFUMX > NCOEFU */
4269 /* or NDIMAX = NDIMEN, NCFUMX = NCOEFU and NCFVMX > NCOEFV */
4271 /* These conditions are not tested in the program. */
4274 /* **********************************************************************
4278 /* Parameter adjustments */
4279 tabini_dim1 = *ndimax;
4280 tabini_dim2 = *ncfumx;
4281 tabini_offset = tabini_dim1 * (tabini_dim2 + 1) + 1;
4282 tabini -= tabini_offset;
4283 tabres_dim1 = *ndimen;
4284 tabres_dim2 = *ncoefu;
4285 tabres_offset = tabres_dim1 * (tabres_dim2 + 1) + 1;
4286 tabres -= tabres_offset;
4289 if (*ndimen == *ndimax) {
4293 /* ----------------------- Compression NDIMEN<>NDIMAX -------------------
4297 for (k = 1; k <= i__1; ++k) {
4299 for (j = 1; j <= i__2; ++j) {
4301 for (i__ = 1; i__ <= i__3; ++i__) {
4302 tabres[i__ + (j + k * tabres_dim2) * tabres_dim1] = tabini[
4303 i__ + (j + k * tabini_dim2) * tabini_dim1];
4312 /* ----------------------- Compression NDIMEN=NDIMAX --------------------
4316 if (*ncoefu == *ncfumx) {
4319 ilong = (*ndimen << 3) * *ncoefu;
4321 for (k = 1; k <= i__1; ++k) {
4322 AdvApp2Var_SysBase::mcrfill_(&ilong,
4323 &tabini[(k * tabini_dim2 + 1) * tabini_dim1 + 1],
4324 &tabres[(k * tabres_dim2 + 1) * tabres_dim1 + 1]);
4329 /* ----------------- Compression NDIMEN=NDIMAX,NCOEFU=NCFUMX ------------
4333 ilong = (*ndimen << 3) * *ncoefu * *ncoefv;
4334 AdvApp2Var_SysBase::mcrfill_(&ilong,
4335 &tabini[tabini_offset],
4336 &tabres[tabres_offset]);
4339 /* ---------------------------- The end ---------------------------------
4346 //=======================================================================
4347 //function : AdvApp2Var_MathBase::mmfmcar_
4349 //=======================================================================
4350 int AdvApp2Var_MathBase::mmfmcar_(integer *ndimen,
4364 /* System generated locals */
4365 integer patold_dim1, patold_dim2, patnew_dim1, patnew_dim2,
4366 i__1, patold_offset,patnew_offset;
4368 /* Local variables */
4369 doublereal* tbaux = 0;
4370 integer ksize, numax, kk;
4374 /* ***********************************************************************
4379 /* LIMITATION OF A SQUARE DEFINED ON (0,1)*(0,1) BETWEEN ISOS */
4380 /* UPARA1 AND UPARA2 (BY U) AND VPARA1 AND VPARA2 BY V. */
4384 /* LIMITATION , SQUARE , PARAMETER */
4386 /* INPUT ARGUMENTS : */
4387 /* ------------------ */
4388 /* NCOFMX: MAX NUMBER OF COEFF OF THE SQUARE BY U */
4389 /* NCOEFU: NUMBER OF COEFF OF THE SQUARE BY U */
4390 /* NCOEFV: NUMBER OF COEFF OF THE SQUARE BY V */
4391 /* PATOLD : THE SQUARE IS LIMITED BY UPARA1,UPARA2 AND VPARA1,VPARA2
4393 /* UPARA1 : LOWER LIMIT OF U */
4394 /* UPARA2 : UPPER LIMIT OF U */
4395 /* VPARA1 : LOWER LIMIT OF V */
4396 /* VPARA2 : UPPER LIMIT OF V */
4398 /* OUTPUT ARGUMENTS : */
4399 /* ------------------- */
4400 /* PATNEW : RELIMITED SQUARE, DEFINED ON (0,1)**2 */
4401 /* IERCOD : =10 COEFF NB TOO GREAT OR NULL */
4402 /* =13 PB IN THE DYNAMIC ALLOCATION */
4405 /* COMMONS USED : */
4406 /* ---------------- */
4408 /* DESCRIPTION/NOTES/LIMITATIONS : */
4409 /* ----------------------------------- */
4410 /* ---> The following call : */
4411 /* CALL MMFMCAR(NCOFMX,NCOEFU,NCOEFV,PATOLD,UPARA1,UPARA2,VPARA1,VPARA2
4414 /* where PATOLD is input/output argument is absolutely legal. */
4416 /* ---> The max number of coeff by u and v of PATOLD is 61 */
4418 /* ---> If NCOEFU < NCOFMX, the data is compressed by MMFMCA9 before
4419 /* limitation by v to get time during the execution */
4420 /* of MMARC41 that follows (the square is processed as a curve of
4422 /* dimension NDIMEN*NCOEFU possessing NCOEFV coefficients). */
4424 /* ***********************************************************************
4427 /* Name of the routine */
4430 /* Parameter adjustments */
4431 patnew_dim1 = *ndimen;
4432 patnew_dim2 = *ncofmx;
4433 patnew_offset = patnew_dim1 * (patnew_dim2 + 1) + 1;
4434 patnew -= patnew_offset;
4435 patold_dim1 = *ndimen;
4436 patold_dim2 = *ncofmx;
4437 patold_offset = patold_dim1 * (patold_dim2 + 1) + 1;
4438 patold -= patold_offset;
4441 ibb = AdvApp2Var_SysBase::mnfndeb_();
4443 AdvApp2Var_SysBase::mgenmsg_("MMFMCAR", 7L);
4447 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
4449 /* **********************************************************************
4451 /* TEST OF COEFFICIENT NUMBERS */
4452 /* **********************************************************************
4455 if (*ncofmx < *ncoefu) {
4459 if (*ncoefu < 1 || *ncoefu > 61 || *ncoefv < 1 || *ncoefv > 61) {
4464 /* **********************************************************************
4466 /* CASE WHEN UPARA1=VPARA1=0 AND UPARA2=VPARA2=1 */
4467 /* **********************************************************************
4470 if (*upara1 == 0. && *upara2 == 1. && *vpara1 == 0. && *vpara2 == 1.) {
4471 ksize = (*ndimen << 3) * *ncofmx * *ncoefv;
4472 AdvApp2Var_SysBase::mcrfill_(&ksize,
4473 &patold[patold_offset],
4474 &patnew[patnew_offset]);
4478 /* **********************************************************************
4480 /* LIMITATION BY U */
4481 /* **********************************************************************
4484 if (*upara1 == 0. && *upara2 == 1.) {
4488 for (kk = 1; kk <= i__1; ++kk) {
4489 mmarc41_(ndimen, ndimen, ncoefu, &patold[(kk * patold_dim2 + 1) *
4490 patold_dim1 + 1], upara1, upara2, &patnew[(kk * patnew_dim2 +
4491 1) * patnew_dim1 + 1], iercod);
4495 /* **********************************************************************
4497 /* LIMITATION BY V */
4498 /* **********************************************************************
4502 if (*vpara1 == 0. && *vpara2 == 1.) {
4506 /* ----------- LIMITATION BY V (WITH COMPRESSION I.E. NCOEFU<NCOFMX) ----
4509 numax = *ndimen * *ncoefu;
4510 if (*ncofmx != *ncoefu) {
4511 /* ------------------------- Dynamic allocation -------------------
4513 ksize = *ndimen * *ncoefu * *ncoefv;
4514 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &ksize, tbaux, &iofst, &ier);
4519 /* --------------- Compression by (NDIMEN,NCOEFU,NCOEFV) ------------
4521 if (*upara1 == 0. && *upara2 == 1.) {
4522 AdvApp2Var_MathBase::mmfmca9_(ndimen,
4528 &patold[patold_offset],
4531 AdvApp2Var_MathBase::mmfmca9_(ndimen,
4537 &patnew[patnew_offset],
4540 /* ------------------------- Limitation by v ------------------------
4542 mmarc41_(&numax, &numax, ncoefv, &tbaux[iofst], vpara1, vpara2, &
4543 tbaux[iofst], iercod);
4544 /* --------------------- Expansion of TBAUX into PATNEW -------------
4546 AdvApp2Var_MathBase::mmfmca8_(ndimen, ncoefu, ncoefv, ndimen, ncofmx, ncoefv, &tbaux[iofst]
4547 , &patnew[patnew_offset]);
4550 /* -------- LIMITATION BY V (WITHOUT COMPRESSION I.E. NCOEFU=NCOFMX) ---
4554 if (*upara1 == 0. && *upara2 == 1.) {
4555 mmarc41_(&numax, &numax, ncoefv, &patold[patold_offset], vpara1,
4556 vpara2, &patnew[patnew_offset], iercod);
4558 mmarc41_(&numax, &numax, ncoefv, &patnew[patnew_offset], vpara1,
4559 vpara2, &patnew[patnew_offset], iercod);
4564 /* **********************************************************************
4567 /* **********************************************************************
4572 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &ksize, tbaux, &iofst, &ier);
4578 /* ------------------------------ The end -------------------------------
4583 AdvApp2Var_SysBase::maermsg_("MMFMCAR", iercod, 7L);
4586 AdvApp2Var_SysBase::mgsomsg_("MMFMCAR", 7L);
4592 //=======================================================================
4593 //function : AdvApp2Var_MathBase::mmfmcb5_
4595 //=======================================================================
4596 int AdvApp2Var_MathBase::mmfmcb5_(integer *isenmsc,
4607 /* System generated locals */
4608 integer courb1_dim1, courb1_offset, courb2_dim1, courb2_offset, i__1,
4611 /* Local variables */
4612 integer i__, nboct, nd;
4615 /* **********************************************************************
4620 /* Reformating (and eventual compression/decompression) of curve */
4621 /* (ndim,.) by (.,ndim) and vice versa. */
4625 /* ALL , MATH_ACCES :: */
4626 /* COURBE&, REORGANISATION,COMPRESSION,INVERSION , &COURBE */
4628 /* INPUT ARGUMENTS : */
4629 /* -------------------- */
4630 /* ISENMSC : required direction of the transfer : */
4631 /* 1 : passage of (NDIMEN,.) ---> (.,NDIMEN) direction to AB
4633 /* -1 : passage of (.,NDIMEN) ---> (NDIMEN,.) direction to TS,T
4635 /* NDIMAX : format / dimension */
4636 /* NCF1MX : format by t of COURB1 */
4637 /* if ISENMSC= 1 : COURB1: The curve to be processed (NDIMAX,.) */
4638 /* NCOEFF : number of coeff of the curve */
4639 /* NCF2MX : format by t of COURB2 */
4640 /* NDIMEN : dimension of the curve and format of COURB2 */
4641 /* if ISENMSC=-1 : COURB2: The curve to be processed (.,NDIMEN) */
4643 /* OUTPUT ARGUMENTS : */
4644 /* --------------------- */
4645 /* if ISENMSC= 1 : COURB2: The resulting curve (.,NDIMEN) */
4646 /* if ISENMSC=-1 : COURB1: The resulting curve (NDIMAX,.) */
4648 /* COMMONS USED : */
4649 /* ------------------ */
4651 /* REFERENCES CALLED : */
4652 /* --------------------- */
4654 /* DESCRIPTION/NOTES/LIMITATIONS : */
4655 /* ----------------------------------- */
4656 /* allow to process the usual transfers as follows : */
4657 /* | ---- ISENMSC = 1 ---- | | ---- ISENMSC =-1 ----- | */
4658 /* TS (3,21) --> (21,3) AB ; AB (21,3) --> (3,21) TS */
4659 /* TS (3,21) --> (NU,3) AB ; AB (NU,3) --> (3,21) TS */
4660 /* (3,NU) --> (21,3) AB ; AB (21,3) --> (3,NU) */
4661 /* (3,NU) --> (NU,3) AB ; AB (NU,3) --> (3,NU) */
4663 /* ***********************************************************************
4667 /* Parameter adjustments */
4668 courb1_dim1 = *ndimax;
4669 courb1_offset = courb1_dim1 + 1;
4670 courb1 -= courb1_offset;
4671 courb2_dim1 = *ncf2mx;
4672 courb2_offset = courb2_dim1 + 1;
4673 courb2 -= courb2_offset;
4676 if (*ndimen > *ndimax || *ncoeff > *ncf1mx || *ncoeff > *ncf2mx) {
4680 if (*ndimen == 1 && *ncf1mx == *ncf2mx) {
4681 nboct = *ncf2mx << 3;
4682 if (*isenmsc == 1) {
4683 AdvApp2Var_SysBase::mcrfill_(&nboct,
4684 &courb1[courb1_offset],
4685 &courb2[courb2_offset]);
4687 if (*isenmsc == -1) {
4688 AdvApp2Var_SysBase::mcrfill_(&nboct,
4689 &courb2[courb2_offset],
4690 &courb1[courb1_offset]);
4697 if (*isenmsc == 1) {
4699 for (nd = 1; nd <= i__1; ++nd) {
4701 for (i__ = 1; i__ <= i__2; ++i__) {
4702 courb2[i__ + nd * courb2_dim1] = courb1[nd + i__ *
4708 } else if (*isenmsc == -1) {
4710 for (nd = 1; nd <= i__1; ++nd) {
4712 for (i__ = 1; i__ <= i__2; ++i__) {
4713 courb1[nd + i__ * courb1_dim1] = courb2[i__ + nd *
4725 /* ***********************************************************************
4733 AdvApp2Var_SysBase::maermsg_("MMFMCB5", iercod, 7L);
4738 //=======================================================================
4739 //function : AdvApp2Var_MathBase::mmfmtb1_
4741 //=======================================================================
4742 int AdvApp2Var_MathBase::mmfmtb1_(integer *maxsz1,
4754 /* System generated locals */
4755 integer table1_dim1, table1_offset, table2_dim1, table2_offset, i__1,
4758 /* Local variables */
4759 doublereal* work = 0;
4760 integer ilong, isize, ii, jj, ier;
4761 intptr_t iofst,iipt, jjpt;
4764 /************************************************************************
4769 /* Inversion of elements of a rectangular table (T1(i,j) */
4770 /* loaded in T2(j,i)) */
4774 /* ALL, MATH_ACCES :: TABLEAU&, INVERSION, &TABLEAU */
4776 /* INPUT ARGUMENTS : */
4777 /* ------------------ */
4778 /* MAXSZ1: Max Nb of elements by the 1st dimension of TABLE1. */
4779 /* TABLE1: Table of reals by two dimensions. */
4780 /* ISIZE1: Nb of useful elements of TABLE1 on the 1st dimension */
4781 /* JSIZE1: Nb of useful elements of TABLE1 on the 2nd dimension */
4782 /* MAXSZ2: Nb max of elements by the 1st dimension of TABLE2. */
4784 /* OUTPUT ARGUMENTS : */
4785 /* ------------------- */
4786 /* TABLE2: Table of reals by two dimensions, containing the transposition
4787 /* of the rectangular table TABLE1. */
4788 /* ISIZE2: Nb of useful elements of TABLE2 on the 1st dimension */
4789 /* JSIZE2: Nb of useful elements of TABLE2 on the 2nd dimension */
4790 /* IERCOD: Erroe coder. */
4792 /* = 1, error in the dimension of tables */
4793 /* ether MAXSZ1 < ISIZE1 (table TABLE1 too small). */
4794 /* or MAXSZ2 < JSIZE1 (table TABLE2 too small). */
4796 /* COMMONS USED : */
4797 /* ---------------- */
4799 /* REFERENCES CALLED : */
4800 /* ---------------------- */
4802 /* DESCRIPTION/NOTES/LIMITATIONS : */
4803 /* ----------------------------------- */
4804 /* It is possible to use TABLE1 as input and output table i.e. */
4806 /* CALL MMFMTB1(MAXSZ1,TABLE1,ISIZE1,JSIZE1,MAXSZ2,TABLE1 */
4807 /* ,ISIZE2,JSIZE2,IERCOD) */
4810 /* **********************************************************************
4814 /* Parameter adjustments */
4815 table1_dim1 = *maxsz1;
4816 table1_offset = table1_dim1 + 1;
4817 table1 -= table1_offset;
4818 table2_dim1 = *maxsz2;
4819 table2_offset = table2_dim1 + 1;
4820 table2 -= table2_offset;
4821 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
4825 if (*isize1 > *maxsz1 || *jsize1 > *maxsz2) {
4830 isize = *maxsz2 * *isize1;
4831 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &isize, work, &iofst, &ier);
4836 /* DO NOT BE AFRAID OF CRUSHING. */
4839 for (ii = 1; ii <= i__1; ++ii) {
4840 iipt = (ii - 1) * *maxsz2 + iofst;
4842 for (jj = 1; jj <= i__2; ++jj) {
4843 jjpt = iipt + (jj - 1);
4844 work[jjpt] = table1[ii + jj * table1_dim1];
4850 AdvApp2Var_SysBase::mcrfill_(&ilong,
4852 &table2[table2_offset]);
4854 /* -------------- The number of elements of TABLE2 is returned ------------
4863 /* ------------------------------- THE END ------------------------------
4865 /* --> Invalid input. */
4869 /* --> Pb of allocation. */
4876 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &isize, work, &iofst, &ier);
4884 //=======================================================================
4885 //function : AdvApp2Var_MathBase::mmgaus1_
4887 //=======================================================================
4888 int AdvApp2Var_MathBase::mmgaus1_(integer *ndimf,
4905 /* System generated locals */
4908 /* Local variables */
4912 doublereal t, u[20], x;
4914 doublereal c1x, c2x;
4915 /* **********************************************************************
4921 /* Calculate the integral of function BFUNX passed in parameter */
4922 /* between limits XD and XF . */
4923 /* The function should be calculated for any value */
4924 /* of the variable in the given interval.. */
4925 /* The method GAUSS-LEGENDRE is used.
4926 /* For explications refer to the book : */
4927 /* Complements de mathematiques a l'usage des Ingenieurs de */
4928 /* l'electrotechnique et des telecommunications. */
4929 /* Par Andre ANGOT - Collection technique et scientifique du CNET
4932 /* The degree of LEGENDRE polynoms used is passed in parameter.
4936 /* INTEGRATION,LEGENDRE,GAUSS */
4938 /* INPUT ARGUMENTS : */
4939 /* ------------------ */
4941 /* NDIMF : Dimension of the function */
4942 /* BFUNX : Function to integrate passed as argument */
4943 /* Should be declared as EXTERNAL in the call routine. */
4944 /* SUBROUTINE BFUNX(NDIMF,X,VAL,IER) */
4946 /* K : Parameter determining the degree of the LEGENDRE polynom that
4948 /* can take a value between 0 and 10. */
4949 /* The degree of the polynom is equal to 4 k, that is 4, 8,
4951 /* 12, 16, 20, 24, 28, 32, 36 and 40. */
4952 /* If K is not correct, the degree is set to 40 directly.
4954 /* XD : Lower limit of the interval of integration. */
4955 /* XF : Upper limit of the interval of integration. */
4956 /* SAUX1 : Auxiliary table */
4957 /* SAUX2 : Auxiliary table */
4959 /* OUTPUT ARGUMENTS : */
4960 /* ------------------- */
4962 /* SOMME : Value of the integral */
4963 /* NITER : Number of iterations to be carried out. */
4964 /* It is equal to the degree of the polynom. */
4966 /* IER : Error code : */
4967 /* < 0 ==> Attention - Warning */
4968 /* = 0 ==> Everything is OK */
4969 /* > 0 ==> Critical error - Apply special processing */
4970 /* ==> Error in the calculation of BFUNX (return code */
4971 /* of this routine */
4973 /* If error => SUM = 0 */
4975 /* COMMONS USED : */
4976 /* ----------------- */
4980 /* REFERENCES CALLED : */
4981 /* ---------------------- */
4984 /* @ BFUNX MVGAUS0 */
4986 /* DESCRIPTION/NOTES/LIMITATIONS : */
4987 /* --------------------------------- */
4989 /* See the explanations detailed in the listing */
4990 /* Use of the GAUSS method (orthogonal polynoms) */
4991 /* The symmetry of roots of these polynomes is used */
4992 /* Depending on K, the degree of the interpolated polynom grows.
4994 /* If you wish to calculate the integral with a given precision, */
4995 /* loop on k varying from 1 to 10 and test the difference of 2
4997 /* consecutive iterations. Stop the loop if this difference is less that
4998 /* an epsilon value set to 10E-6 for example. */
4999 /* If S1 and S2 are 2 successive iterations, test following this example :
5002 /* AF=DABS(S1-S2) */
5004 /* If AS < 1 test if FS < eps otherwise test if AF/AS < eps
5006 /* -- ----- ----- */
5008 /************************************************************************
5011 /************************************************************************
5016 /* ****** General Initialization */
5018 /* Parameter adjustments */
5024 AdvApp2Var_SysBase::mvriraz_(ndimf,
5028 /* ****** Loading of coefficients U and H ** */
5029 /* -------------------------------------------- */
5031 mvgaus0_(k, u, h__, &ndeg, iercod);
5036 /* ****** C1X => Medium interval point [XD,XF] */
5037 /* ****** C2X => 1/2 amplitude interval [XD,XF] */
5039 c1x = (*xf + *xd) * .5;
5040 c2x = (*xf - *xd) * .5;
5042 /* ---------------------------------------- */
5043 /* ****** Integration for degree NDEG ** */
5044 /* ---------------------------------------- */
5047 for (j = 1; j <= i__1; ++j) {
5051 (*bfunx)(ndimf, &x, &saux1[1], iercod);
5057 (*bfunx)(ndimf, &x, &saux2[1], iercod);
5063 for (idimf = 1; idimf <= i__2; ++idimf) {
5064 somme[idimf] += h__[j - 1] * (saux1[idimf] + saux2[idimf]);
5071 for (idimf = 1; idimf <= i__1; ++idimf) {
5072 somme[idimf] *= c2x;
5075 /* ****** End of sub-program ** */
5081 //=======================================================================
5082 //function : mmherm0_
5084 //=======================================================================
5085 int mmherm0_(doublereal *debfin,
5088 integer c__576 = 576;
5092 /* System generated locals */
5096 /* Local variables */
5097 doublereal amat[36] /* was [6][6] */;
5100 integer iord1, iord2;
5101 doublereal miden[36] /* was [6][6] */;
5103 doublereal epspi, d1, d2;
5104 integer ii, jj, pp, ncf;
5106 integer iof[2], ier;
5107 doublereal mat[36] /* was [6][6] */;
5109 doublereal abid[72] /* was [12][6] */;
5110 /* ***********************************************************************
5115 /* INIT OF COEFFS. OF POLYNOMS OF HERMIT INTERPOLATION */
5119 /* MATH_ACCES :: HERMITE */
5121 /* INPUT ARGUMENTS */
5122 /* -------------------- */
5123 /* DEBFIN : PARAMETERS DEFINING THE CONSTRAINTS */
5124 /* DEBFIN(1) : FIRST PARAMETER */
5125 /* DEBFIN(2) : SECOND PARAMETER */
5127 /* ONE SHOULD HAVE: */
5128 /* ABS (DEBFIN(I)) < 100 */
5130 /* (ABS(DEBFIN(1)+ABS(DEBFIN(2))) > 1/100 */
5131 /* (for overflows) */
5133 /* ABS(DEBFIN(2)-DEBFIN(1)) / (ABS(DEBFIN(1)+ABS(DEBFIN(2))) > 1/100
5135 /* (for the conditioning) */
5138 /* OUTPUT ARGUMENTS : */
5139 /* --------------------- */
5141 /* IERCOD : Error code : 0 : O.K. */
5142 /* 1 : value of DEBFIN */
5143 /* are unreasonable */
5144 /* -1 : init was already done */
5145 /* (OK but no processing) */
5147 /* COMMONS USED : */
5148 /* ------------------ */
5150 /* REFERENCES CALLED : */
5151 /* ---------------------- */
5154 /* DESCRIPTION/NOTES/LIMITATIONS : */
5155 /* ----------------------------------- */
5157 /* This program initializes the coefficients of Hermit polynoms */
5158 /* that are read later by MMHERM1 */
5159 /* ***********************************************************************
5164 /* **********************************************************************
5169 /* Used to STORE coefficients of Hermit interpolation polynoms
5175 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5176 /* ----------------------------------- */
5178 /* The coefficients of hermit polynoms are calculated by */
5179 /* the routine MMHERM0 and read by the routine MMHERM1 */
5181 /* **********************************************************************
5188 /* NBCOEF is the size of CMHERM (see below) */
5189 /* ***********************************************************************
5198 /* ***********************************************************************
5201 /* ***********************************************************************
5205 /* Parameter adjustments */
5209 d1 = advapp_abs(debfin[1]);
5210 if (d1 > (float)100.) {
5214 d2 = advapp_abs(debfin[2]);
5215 if (d2 > (float)100.) {
5220 if (d2 < (float).01) {
5224 d1 = (d__1 = debfin[2] - debfin[1], advapp_abs(d__1));
5225 if (d1 / d2 < (float).01) {
5230 /* ***********************************************************************
5232 /* Initialization */
5233 /* ***********************************************************************
5241 /* ***********************************************************************
5244 /* IS IT ALREADY INITIALIZED ? */
5246 d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
5249 if (debfin[1] != mmcmher_.tdebut) {
5252 if (debfin[2] != mmcmher_.tfinal) {
5255 if (d1 != mmcmher_.verifi) {
5263 /* ***********************************************************************
5266 /* ***********************************************************************
5272 /* Init. matrix identity : */
5275 AdvApp2Var_SysBase::mvriraz_(&ncmat,
5278 for (ii = 1; ii <= 6; ++ii) {
5279 miden[ii + ii * 6 - 7] = 1.;
5285 /* Init to 0 of table CMHERM */
5287 AdvApp2Var_SysBase::mvriraz_(&c__576, mmcmher_.cmherm);
5289 /* Calculation by solution of linear systems */
5291 for (iord1 = -1; iord1 <= 2; ++iord1) {
5292 for (iord2 = -1; iord2 <= 2; ++iord2) {
5299 iof[1] = iord[0] + 1;
5302 ncf = iord[0] + iord[1] + 2;
5304 /* Calculate matrix MAT to invert: */
5306 for (cot = 1; cot <= 2; ++cot) {
5309 if (iord[cot - 1] > -1) {
5312 for (jj = 1; jj <= i__1; ++jj) {
5318 i__1 = iord[cot - 1] + 1;
5319 for (pp = 1; pp <= i__1; ++pp) {
5321 ii = pp + iof[cot - 1];
5326 for (jj = 1; jj <= i__2; ++jj) {
5327 mat[ii + jj * 6 - 7] = (float)0.;
5332 for (jj = pp; jj <= i__2; ++jj) {
5334 /* everything is done in these 3 lines
5337 mat[ii + jj * 6 - 7] = cof[jj - 1] * prod;
5338 cof[jj - 1] *= jj - pp;
5339 prod *= debfin[cot];
5352 AdvApp2Var_MathBase::mmmrslwd_(&c__6, &ncf, &ncf, mat, miden, &epspi, abid, amat, &
5359 for (cot = 1; cot <= 2; ++cot) {
5360 i__1 = iord[cot - 1] + 1;
5361 for (pp = 1; pp <= i__1; ++pp) {
5363 for (ii = 1; ii <= i__2; ++ii) {
5364 mmcmher_.cmherm[ii + (pp + (cot + ((iord1 + (iord2 <<
5365 2)) << 1)) * 3) * 6 + 155] = amat[ii + (pp +
5366 iof[cot - 1]) * 6 - 7];
5379 /* ***********************************************************************
5382 /* The initialized flag is located: */
5384 mmcmher_.tdebut = debfin[1];
5385 mmcmher_.tfinal = debfin[2];
5387 d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
5388 mmcmher_.verifi = d1 * 16111959;
5391 /* ***********************************************************************
5396 /* ***********************************************************************
5407 /* ***********************************************************************
5412 AdvApp2Var_SysBase::maermsg_("MMHERM0", iercod, 7L);
5414 /* ***********************************************************************
5419 //=======================================================================
5420 //function : mmherm1_
5422 //=======================================================================
5423 int mmherm1_(doublereal *debfin,
5429 /* System generated locals */
5430 integer hermit_dim1, hermit_dim2, hermit_offset;
5432 /* Local variables */
5437 /* ***********************************************************************
5442 /* reading of coeffs. of HERMIT interpolation polynoms */
5446 /* MATH_ACCES :: HERMIT */
5448 /* INPUT ARGUMENTS : */
5449 /* -------------------- */
5450 /* DEBFIN : PARAMETES DEFINING THE CONSTRAINTS */
5451 /* DEBFIN(1) : FIRST PARAMETER */
5452 /* DEBFIN(2) : SECOND PARAMETER */
5454 /* Should be equal to the corresponding arguments during the */
5455 /* last call to MMHERM0 for the initialization of coeffs. */
5457 /* ORDRMX : indicates the dimensioning of HERMIT: */
5458 /* there is no choice : ORDRMX should be equal to the value */
5459 /* of PARAMETER IORDMX of INCLUDE MMCMHER, or 2 for the moment */
5461 /* IORDRE (2) : Orders of constraints in each corresponding parameter DEBFIN(I)
5462 /* should be between -1 (no constraints) and ORDRMX. */
5465 /* OUTPUT ARGUMENTS : */
5466 /* --------------------- */
5468 /* HERMIT : HERMIT(1:IORDRE(1)+IORDRE(2)+2, j, cote) are the */
5469 /* coefficients in the canonic base of Hermit polynom */
5470 /* corresponding to orders IORDRE with parameters DEBFIN for */
5471 /* the constraint of order j on DEBFIN(cote). j is between 0 and IORDRE(cote). */
5474 /* IERCOD : Error code : */
5475 /* -1: O.K but necessary to reinitialize the coefficients */
5476 /* (info for optimization) */
5478 /* 1 : Error in MMHERM0 */
5479 /* 2 : arguments invalid */
5481 /* COMMONS USED : */
5482 /* ------------------ */
5484 /* REFERENCES CALLED : */
5485 /* ---------------------- */
5488 /* DESCRIPTION/NOTES/LIMITATIONS : */
5489 /* ----------------------------------- */
5491 /* This program reads coefficients of Hermit polynoms */
5492 /* that were earlier initialized by MMHERM0 */
5494 /* PMN : initialisation is no more done by the caller. */
5497 /* ***********************************************************************
5502 /* **********************************************************************
5507 /* Serves to STORE the coefficients of Hermit interpolation polynoms
5513 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5514 /* ----------------------------------- */
5516 /* the coefficients of Hetmit polynoms are calculated by */
5517 /* routine MMHERM0 and read by routine MMHERM1 */
5520 /* **********************************************************************
5527 /* NBCOEF is the size of CMHERM (see lower) */
5531 /* ***********************************************************************
5538 /* ***********************************************************************
5540 /* Initializations */
5541 /* ***********************************************************************
5544 /* Parameter adjustments */
5546 hermit_dim1 = (*ordrmx << 1) + 2;
5547 hermit_dim2 = *ordrmx + 1;
5548 hermit_offset = hermit_dim1 * hermit_dim2 + 1;
5549 hermit -= hermit_offset;
5556 /* ***********************************************************************
5559 /* ***********************************************************************
5567 for (cot = 1; cot <= 2; ++cot) {
5568 if (iordre[cot] < -1) {
5571 if (iordre[cot] > *ordrmx) {
5578 /* IS-IT CORRECTLY INITIALIZED ? */
5580 d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
5583 /* OTHERWISE IT IS INITIALIZED */
5585 if (debfin[1] != mmcmher_.tdebut || debfin[2] != mmcmher_.tfinal || d1
5586 != mmcmher_.verifi) {
5588 mmherm0_(&debfin[1], iercod);
5595 /* ***********************************************************************
5598 /* ***********************************************************************
5603 AdvApp2Var_SysBase::msrfill_(&nbval, &mmcmher_.cmherm[((((iordre[1] + (iordre[2] << 2)) << 1)
5604 + 1) * 3 + 1) * 6 + 156], &hermit[hermit_offset]);
5606 /* ***********************************************************************
5611 /* ***********************************************************************
5622 /* ***********************************************************************
5627 AdvApp2Var_SysBase::maermsg_("MMHERM1", iercod, 7L);
5629 /* ***********************************************************************
5634 //=======================================================================
5635 //function : AdvApp2Var_MathBase::mmhjcan_
5637 //=======================================================================
5638 int AdvApp2Var_MathBase::mmhjcan_(integer *ndimen,
5651 /* System generated locals */
5652 integer tcbold_dim1, tcbold_dim2, tcbold_offset, tcbnew_dim1, tcbnew_dim2,
5653 tcbnew_offset, i__1, i__2, i__3, i__4, i__5;
5656 /* Local variables */
5659 doublereal taux1[21];
5660 integer d__, e, i__, k;
5663 doublereal tjacap[21];
5665 doublereal hermit[36]/* was [6][3][2] */, ctenor, bornes[2];
5669 /* ***********************************************************************
5674 /* CONVERSION OF TABLE TCBOLD OF POLYNOMIAL CURVE COEFFICIENTS */
5675 /* EXPRESSED IN HERMIT JACOBI BASE, INTO A */
5676 /* TABLE OF COEFFICIENTS TCBNEW OF COURVES EXPRESSED IN THE CANONIC BASE */
5680 /* CANNONIC, HERMIT, JACCOBI */
5682 /* INPUT ARGUMENTS : */
5683 /* -------------------- */
5684 /* ORDHER : ORDER OF HERMIT POLYNOMS OR ORDER OF CONTINUITY */
5685 /* NCOEFS : NUMBER OF COEFFICIENTS OF A POLYNOMIAL CURVE */
5686 /* FOR ONE OF ITS NDIM COMPONENTS;(DEGREE+1 OF THE CURVE)
5688 /* NDIM : DIMENSION OF THE CURVE */
5689 /* CBHEJA : TABLE OF COEFFICIENTS OF THE CURVE IN THE BASE */
5691 /* (H(0,-1),..,H(ORDHER,-1),H(0,1),..,H(ORDHER,1), */
5692 /* JA(ORDHER+1,2*ORDHER+2),....,JA(ORDHER+1,NCOEFS-1) */
5694 /* OUTPUT ARGUMENTS : */
5695 /* --------------------- */
5696 /* CBRCAN : TABLE OF COEFFICIENTS OF THE CURVE IN THE CANONIC BASE */
5699 /* COMMONS USED : */
5700 /* ------------------ */
5703 /* REFERENCES CALLED : */
5704 /* --------------------- */
5707 /* ***********************************************************************
5711 /* ***********************************************************************
5716 /* Providesinteger constants from 0 to 1000 */
5722 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5723 /* ----------------------------------- */
5725 /* ***********************************************************************
5729 /* ***********************************************************************
5735 /* ***********************************************************************
5737 /* INITIALIZATION */
5738 /* ***********************************************************************
5741 /* Parameter adjustments */
5743 tcbnew_dim1 = *ndimen;
5744 tcbnew_dim2 = *ncflim;
5745 tcbnew_offset = tcbnew_dim1 * (tcbnew_dim2 + 1) + 1;
5746 tcbnew -= tcbnew_offset;
5747 tcbold_dim1 = *ndimen;
5748 tcbold_dim2 = *ncflim;
5749 tcbold_offset = tcbold_dim1 * (tcbold_dim2 + 1) + 1;
5750 tcbold -= tcbold_offset;
5753 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
5755 AdvApp2Var_SysBase::mgenmsg_("MMHJCAN", 7L);
5762 /* ***********************************************************************
5765 /* ***********************************************************************
5775 /* CALCULATION OF HERMIT POLYNOMS IN THE CANONIC BASE ON (-1,1) */
5778 iordre[0] = *orcont;
5779 iordre[1] = *orcont;
5780 mmherm1_(bornes, &c__2, iordre, hermit, &ier);
5790 for (e = 1; e <= i__1; ++e) {
5792 ctenor = (tdecop[e] - tdecop[e - 1]) / 2;
5800 for (d__ = 1; d__ <= i__2; ++d__) {
5802 /* CONVERSION OF THE COEFFICIENTS OF THE PART OF THE CURVE EXPRESSED */
5803 /* IN HERMIT BASE, INTO THE CANONIC BASE */
5805 AdvApp2Var_SysBase::mvriraz_(&ncoeff, taux1);
5808 for (k = 1; k <= i__3; ++k) {
5810 for (i__ = 1; i__ <= i__4; ++i__) {
5812 mfact = AdvApp2Var_MathBase::pow__di(&ctenor, &i__5);
5813 taux1[k - 1] += (tcbold[d__ + (i__ + e * tcbold_dim2) *
5814 tcbold_dim1] * hermit[k + (i__ + 2) * 6 - 19] +
5815 tcbold[d__ + (i__ + aux1 + e * tcbold_dim2) *
5816 tcbold_dim1] * hermit[k + (i__ + 5) * 6 - 19]) *
5823 for (i__ = aux2 + 1; i__ <= i__3; ++i__) {
5824 taux1[i__ - 1] = tcbold[d__ + (i__ + e * tcbold_dim2) *
5828 /* CONVERSION OF THE COEFFICIENTS OF THE PART OF THE CURVE EXPRESSED */
5829 /* IN CANONIC-JACOBI BASE, INTO THE CANONIC BASE */
5833 AdvApp2Var_MathBase::mmapcmp_(&minombr_.nbr[1], &c__21, &ncoeff, taux1, tjacap);
5834 AdvApp2Var_MathBase::mmjacan_(orcont, &ndeg, tjacap, taux1);
5836 /* RECOPY THE COEFS RESULTING FROM THE CONVERSION IN THE TABLE */
5840 for (i__ = 1; i__ <= i__3; ++i__) {
5841 tcbnew[d__ + (i__ + e * tcbnew_dim2) * tcbnew_dim1] = taux1[
5850 /* ***********************************************************************
5852 /* PROCESSING OF ERRORS */
5853 /* ***********************************************************************
5863 /* ***********************************************************************
5865 /* RETURN CALLING PROGRAM */
5866 /* ***********************************************************************
5871 AdvApp2Var_SysBase::maermsg_("MMHJCAN", iercod, 7L);
5873 AdvApp2Var_SysBase::mgsomsg_("MMHJCAN", 7L);
5878 //=======================================================================
5879 //function : AdvApp2Var_MathBase::mminltt_
5881 //=======================================================================
5882 int AdvApp2Var_MathBase::mminltt_(integer *ncolmx,
5888 doublereal *,//epseg,
5891 /* System generated locals */
5892 integer tabtri_dim1, tabtri_offset, i__1, i__2;
5894 /* Local variables */
5896 integer icol, ilgn, nlgn, noct, inser;
5897 doublereal epsega = 0.;
5900 /* ***********************************************************************
5905 /* . Insert a line in a table parsed without redundance */
5909 /* TOUS,MATH_ACCES :: TABLEAU&,INSERTION,&TABLEAU */
5911 /* INPUT ARGUMENTS : */
5912 /* -------------------- */
5913 /* . NCOLMX : Number of columns in the table */
5914 /* . NLGNMX : Number of lines in the table */
5915 /* . TABTRI : Table parsed by lines without redundances */
5916 /* . NBRCOL : Number of columns used */
5917 /* . NBRLGN : Number of lines used */
5918 /* . AJOUTE : Line to be added */
5919 /* . EPSEGA : Epsilon to test the redundance */
5921 /* OUTPUT ARGUMENTS : */
5922 /* --------------------- */
5923 /* . TABTRI : Table parsed by lines without redundances */
5924 /* . NBRLGN : Number of lines used */
5925 /* . IERCOD : 0 -> No problem */
5926 /* 1 -> The table is full */
5928 /* COMMONS USED : */
5929 /* ------------------ */
5931 /* REFERENCES CALLED : */
5932 /* --------------------- */
5934 /* DESCRIPTION/NOTES/LIMITATIONS : */
5935 /* ----------------------------------- */
5936 /* . The line is inserted only if there is no line with all
5938 /* elements equl to those which are planned to be insered, to epsilon. */
5940 /* . Level of de debug = 3 */
5944 /* DECLARATIONS , CONTROL OF INPUT ARGUMENTS , INITIALIZATION */
5945 /* ***********************************************************************
5948 /* --- Parameters */
5954 /* --- Local variables */
5959 /* Parameter adjustments */
5960 tabtri_dim1 = *ncolmx;
5961 tabtri_offset = tabtri_dim1 + 1;
5962 tabtri -= tabtri_offset;
5966 ibb = AdvApp2Var_SysBase::mnfndeb_();
5969 AdvApp2Var_SysBase::mgenmsg_("MMINLTT", 7L);
5972 /* --- Control arguments */
5974 if (*nbrlgn >= *nlgnmx) {
5978 /* -------------------- */
5979 /* *** INITIALIZATION */
5980 /* -------------------- */
5984 /* ---------------------------- */
5985 /* *** SEARCH OF REDUNDANCE */
5986 /* ---------------------------- */
5989 for (ilgn = 1; ilgn <= i__1; ++ilgn) {
5990 if (tabtri[ilgn * tabtri_dim1 + 1] >= ajoute[1] - epsega) {
5991 if (tabtri[ilgn * tabtri_dim1 + 1] <= ajoute[1] + epsega) {
5993 for (icol = 1; icol <= i__2; ++icol) {
5994 if (tabtri[icol + ilgn * tabtri_dim1] < ajoute[icol] -
5995 epsega || tabtri[icol + ilgn * tabtri_dim1] >
5996 ajoute[icol] + epsega) {
6010 /* ----------------------------------- */
6011 /* *** SEARCH OF THE INSERTION POINT */
6012 /* ----------------------------------- */
6017 for (ilgn = 1; ilgn <= i__1; ++ilgn) {
6019 for (icol = 1; icol <= i__2; ++icol) {
6020 if (tabtri[icol + ilgn * tabtri_dim1] < ajoute[icol]) {
6023 if (tabtri[icol + ilgn * tabtri_dim1] > ajoute[icol]) {
6034 /* -------------- */
6036 /* -------------- */
6043 /* --- Shift lower */
6045 nlgn = *nbrlgn - inser;
6047 noct = (*ncolmx << 3) * nlgn;
6048 AdvApp2Var_SysBase::mcrfill_(&noct,
6049 &tabtri[inser * tabtri_dim1 + 1],
6050 &tabtri[(inser + 1)* tabtri_dim1 + 1]);
6055 noct = *nbrcol << 3;
6056 AdvApp2Var_SysBase::mcrfill_(&noct,
6058 &tabtri[inser * tabtri_dim1 + 1]);
6062 /* ******************************************************************** */
6063 /* OUTPUT ERROR , RETURN CALLING PROGRAM , MESSAGES */
6064 /* ******************************************************************** */
6066 /* --- The table is already full */
6075 AdvApp2Var_SysBase::maermsg_("MMINLTT", iercod, 7L);
6078 AdvApp2Var_SysBase::mgsomsg_("MMINLTT", 7L);
6083 //=======================================================================
6084 //function : AdvApp2Var_MathBase::mmjacan_
6086 //=======================================================================
6087 int AdvApp2Var_MathBase::mmjacan_(const integer *ideriv,
6092 /* System generated locals */
6093 integer poljac_dim1, i__1, i__2;
6095 /* Local variables */
6096 integer iptt, i__, j, ibb;
6099 /* ***********************************************************************
6104 /* Routine of transfer of Jacobi normalized to canonic [-1,1], */
6105 /* the tables are ranked by even, then by uneven degree. */
6109 /* LEGENDRE,JACOBI,PASSAGE. */
6111 /* INPUT ARGUMENTS : */
6112 /* ------------------ */
6113 /* IDERIV : Order of Jacobi between -1 and 2. */
6114 /* NDEG : The true degree of the polynom. */
6115 /* POLJAC : The polynom in the Jacobi base. */
6117 /* OUTPUT ARGUMENTS : */
6118 /* ------------------- */
6119 /* POLCAN : The curve expressed in the canonic base [-1,1]. */
6121 /* COMMONS USED : */
6122 /* ---------------- */
6124 /* REFERENCES CALLED : */
6125 /* ----------------------- */
6127 /* DESCRIPTION/NOTES/LIMITATIONS : */
6128 /* ----------------------------------- */
6131 /* ***********************************************************************
6134 /* Name of the routine */
6136 /* Matrices of conversion */
6139 /* ***********************************************************************
6144 /* MATRIX OF TRANSFORMATION OF LEGENDRE BASE */
6150 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
6151 /* ----------------------------------- */
6154 /* ***********************************************************************
6159 /* Legendre common / Restricted Casteljau. */
6161 /* 0:1 0 Concerns the even terms, 1 the uneven terms. */
6162 /* CANPLG : Matrix of passage to canonic from Jacobi with calculated parities */
6163 /* PLGCAN : Matrix of passage from Jacobi to canonic with calculated parities */
6166 /* ***********************************************************************
6169 /* Parameter adjustments */
6170 poljac_dim1 = *ndeg / 2 + 1;
6173 ibb = AdvApp2Var_SysBase::mnfndeb_();
6175 AdvApp2Var_SysBase::mgenmsg_("MMJACAN", 7L);
6178 /* ----------------- Expression of terms of even degree ----------------
6182 for (i__ = 0; i__ <= i__1; ++i__) {
6184 iptt = i__ * 31 - (i__ + 1) * i__ / 2 + 1;
6186 for (j = i__; j <= i__2; ++j) {
6187 bid += mmjcobi_.plgcan[iptt + j + *ideriv * 992 + 991] * poljac[
6191 polcan[i__ * 2] = bid;
6195 /* --------------- Expression of terms of uneven degree ----------------
6202 i__1 = (*ndeg - 1) / 2;
6203 for (i__ = 0; i__ <= i__1; ++i__) {
6205 iptt = i__ * 31 - (i__ + 1) * i__ / 2 + 1;
6206 i__2 = (*ndeg - 1) / 2;
6207 for (j = i__; j <= i__2; ++j) {
6208 bid += mmjcobi_.plgcan[iptt + j + ((*ideriv << 1) + 1) * 496 +
6209 991] * poljac[j + poljac_dim1];
6212 polcan[(i__ << 1) + 1] = bid;
6216 /* -------------------------------- The end -----------------------------
6221 AdvApp2Var_SysBase::mgsomsg_("MMJACAN", 7L);
6226 //=======================================================================
6227 //function : AdvApp2Var_MathBase::mmjaccv_
6229 //=======================================================================
6230 int AdvApp2Var_MathBase::mmjaccv_(const integer *ncoef,
6231 const integer *ndim,
6232 const integer *ider,
6233 const doublereal *crvlgd,
6238 /* Initialized data */
6240 static char nomprg[8+1] = "MMJACCV ";
6242 /* System generated locals */
6243 integer crvlgd_dim1, crvlgd_offset, crvcan_dim1, crvcan_offset,
6244 polaux_dim1, i__1, i__2;
6246 /* Local variables */
6247 integer ndeg, i__, nd, ii, ibb;
6249 /* ***********************************************************************
6254 /* Passage from the normalized Jacobi base to the canonic base. */
6258 /* SMOOTHING, BASE, LEGENDRE */
6261 /* INPUT ARGUMENTS : */
6262 /* ------------------ */
6263 /* NDIM: Space Dimension. */
6264 /* NCOEF: Degree +1 of the polynom. */
6265 /* IDER: Order of Jacobi polynoms. */
6266 /* CRVLGD : Curve in the base of Jacobi. */
6268 /* OUTPUT ARGUMENTS : */
6269 /* ------------------- */
6270 /* POLAUX : Auxilliary space. */
6271 /* CRVCAN : The curve in the canonic base [-1,1] */
6273 /* COMMONS USED : */
6274 /* ---------------- */
6276 /* REFERENCES CALLED : */
6277 /* ----------------------- */
6279 /* DESCRIPTION/NOTES/LIMITATIONS : */
6280 /* ----------------------------------- */
6283 /* *********************************************************************
6286 /* Name of the routine */
6287 /* Parameter adjustments */
6288 polaux_dim1 = (*ncoef - 1) / 2 + 1;
6289 crvcan_dim1 = *ncoef - 1 + 1;
6290 crvcan_offset = crvcan_dim1;
6291 crvcan -= crvcan_offset;
6292 crvlgd_dim1 = *ncoef - 1 + 1;
6293 crvlgd_offset = crvlgd_dim1;
6294 crvlgd -= crvlgd_offset;
6298 ibb = AdvApp2Var_SysBase::mnfndeb_();
6300 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
6306 for (nd = 1; nd <= i__1; ++nd) {
6307 /* Loading of the auxilliary table. */
6310 for (i__ = 0; i__ <= i__2; ++i__) {
6311 polaux[i__] = crvlgd[ii + nd * crvlgd_dim1];
6318 i__2 = (ndeg - 1) / 2;
6319 for (i__ = 0; i__ <= i__2; ++i__) {
6320 polaux[i__ + polaux_dim1] = crvlgd[ii + nd * crvlgd_dim1];
6325 /* Call the routine of base change. */
6326 AdvApp2Var_MathBase::mmjacan_(ider, &ndeg, polaux, &crvcan[nd * crvcan_dim1]);
6335 //=======================================================================
6336 //function : mmloncv_
6338 //=======================================================================
6339 int mmloncv_(integer *ndimax,
6349 /* Initialized data */
6353 /* System generated locals */
6354 integer courbe_dim1, courbe_offset, i__1, i__2;
6356 /* Local variables */
6359 doublereal c1, c2, d1, d2, wgaus[20], uroot[20], x1, x2, dd;
6362 doublereal der1, der2;
6367 /* **********************************************************************
6370 /* FUNCTION : Length of an arc of curve on a given interval */
6371 /* ---------- for a function the mathematic representation */
6372 /* which of is a multidimensional polynom. */
6373 /* The polynom is a set of polynoms the coefficients which of are ranked
6374 /* in a table with 2 indices, each line relative to 1 polynom. */
6375 /* The polynom is defined by its coefficients ordered by increasing
6376 * power of the variable. */
6377 /* All polynoms have the same number of coefficients (and the same degree). */
6379 /* KEYWORDS : LENGTH, CURVE */
6382 /* INPUT ARGUMENTS : */
6383 /* -------------------- */
6385 /* NDIMAX : Max number of lines of tables (max number of polynoms). */
6386 /* NDIMEN : Dimension of the polynom (Nomber of polynoms). */
6387 /* NCOEFF : Number of coefficients of the polynom (no limitation) */
6388 /* This is degree + 1 */
6389 /* COURBE : Coefficients of the polynom ordered by increasing power */
6390 /* Dimension to (NDIMAX,NCOEFF). */
6391 /* TDEBUT : Lower limit of integration for length calculation. */
6392 /* TFINAL : Upper limit of integration for length calculation. */
6394 /* OUTPUT ARGUMENTS : */
6395 /* --------------------- */
6396 /* XLONGC : Length of arc of curve */
6398 /* IERCOD : Error code : */
6399 /* = 0 ==> All is OK */
6400 /* = 1 ==> NDIMEN or NCOEFF negative or null */
6401 /* = 2 ==> Pb loading Legendre roots and Gauss weight */
6404 /* If error => XLONGC = 0 */
6406 /* COMMONS USED : */
6407 /* ------------------ */
6411 /* REFERENCES CALLED : */
6412 /* ---------------------- */
6414 /* MAERMSG R*8 DSQRT I*4 MIN */
6417 /* DESCRIPTION/NOTES/LIMITATIONS : */
6418 /* ----------------------------------- */
6420 /* See VGAUSS to understand well the technique. */
6421 /* Actually SQRT (dpi^2) is integrated for i=1,nbdime */
6422 /* Calculation of the derivative is included in the code to avoid an additional */
6423 /* call of the routine. */
6425 /* The integrated function is strictly increasing, it */
6426 /* is not necessary to use a high degree for the GAUSS method GAUSS. */
6428 /* The degree of LEGENDRE polynom results from the degree of the */
6429 /* polynom to be integrated. It can vary from 4 to 40 (with step of 4). */
6431 /* The precision (relative) of integration is of order 1.D-8. */
6433 /* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
6435 /* Attention : the precision of the result is not controlled. */
6436 /* If you wish to control it, use MMCGLC1, taking into account that */
6437 /* the performance (in time) will be worse. */
6439 /* >=====================================================================
6442 /* ATTENTION : SAVE KGAR WGAUS and UROOT EVENTUALLY */
6444 /* INTEGER I1,I20 */
6445 /* PARAMETER (I1=1,I20=20) */
6447 /* Parameter adjustments */
6448 courbe_dim1 = *ndimax;
6449 courbe_offset = courbe_dim1 + 1;
6450 courbe -= courbe_offset;
6454 /* ****** General initialization ** */
6459 /* ****** Initialization of UROOT, WGAUS, NGAUS and KGAR ** */
6461 /* CALL MXVINIT(IERXV,'INTEGER',I1,KGAR,'INTEGER',I1,NGAUS */
6462 /* 1 ,'DOUBLE PRECISION',I20,UROOT,'DOUBLE PRECISION',I20,WGAUS) */
6463 /* IF (IERXV.GT.0) KGAR=0 */
6465 /* ****** Test the equity of limits ** */
6467 if (*tdebut == *tfinal) {
6472 /* ****** Test the dimension and the number of coefficients ** */
6474 if (*ndimen <= 0 || *ncoeff <= 0) {
6479 /* ****** Calculate the optimal degree ** */
6481 kk = *ncoeff / 4 + 1;
6482 kk = advapp_min(kk,10);
6484 /* ****** Return the coefficients for the integral (DEGRE=4*KK) */
6485 /* if KK <> KGAR. */
6488 mvgaus0_(&kk, uroot, wgaus, &ngaus, iercod);
6497 /* C1 => Point medium interval */
6498 /* C2 => 1/2 amplitude interval */
6500 c1 = (*tfinal + *tdebut) * .5;
6501 c2 = (*tfinal - *tdebut) * .5;
6503 /* ----------------------------------------------------------- */
6504 /* ****** Integration - Loop on GAUSS intervals ** */
6505 /* ----------------------------------------------------------- */
6510 for (jj = 1; jj <= i__1; ++jj) {
6512 /* ****** Integration taking the symmetry into account ** */
6514 tran = c2 * uroot[jj - 1];
6518 /* ****** Derivation on the dimension of the space ** */
6523 for (kk = 1; kk <= i__2; ++kk) {
6524 d1 = (*ncoeff - 1) * courbe[kk + *ncoeff * courbe_dim1];
6526 for (ii = *ncoeff - 1; ii >= 2; --ii) {
6527 dd = (ii - 1) * courbe[kk + ii * courbe_dim1];
6537 /* ****** Integration ** */
6539 som += wgaus[jj - 1] * c2 * (sqrt(der1) + sqrt(der2));
6541 /* ****** End of loop on GAUSS intervals ** */
6546 /* ****** Work ended ** */
6550 /* ****** It is forced IERCOD = 0 ** */
6554 /* ****** Final processing ** */
6558 /* ****** Save UROOT, WGAUS, NGAUS and KGAR ** */
6560 /* CALL MXVSAVE(IERXV,'INTEGER',I1,KGAR,'INTEGER',I1,NGAUS */
6561 /* 1 ,'DOUBLE PRECISION',I20,UROOT,'DOUBLE PRECISION',I20,WGAUS) */
6562 /* IF (IERXV.GT.0) KGAR=0 */
6564 /* ****** End of sub-program ** */
6567 AdvApp2Var_SysBase::maermsg_("MMLONCV", iercod, 7L);
6572 //=======================================================================
6573 //function : AdvApp2Var_MathBase::mmpobas_
6575 //=======================================================================
6576 int AdvApp2Var_MathBase::mmpobas_(doublereal *tparam,
6588 /* Initialized data */
6590 doublereal moin11[2] = { -1.,1. };
6592 /* System generated locals */
6593 integer valbas_dim1, i__1;
6595 /* Local variables */
6596 doublereal vjac[80], herm[24];
6599 integer nwcof, iunit;
6600 doublereal wpoly[7];
6601 integer ii, jj, iorjac;
6602 doublereal hermit[36] /* was [6][3][2] */;
6603 integer kk1, kk2, kk3;
6607 /* ***********************************************************************
6612 /* Position on the polynoms of base hermit-Jacobi */
6613 /* and their succesive derivatives */
6617 /* PUBLIC, POSITION, HERMIT, JACOBI */
6619 /* INPUT ARGUMENTS : */
6620 /* -------------------- */
6621 /* TPARAM : Parameter for which the position is found. */
6622 /* IORDRE : Orderof hermit-Jacobi (-1,0,1, ou 2) */
6623 /* NCOEFF : Number of coefficients of polynoms (Nb of value to calculate) */
6624 /* NDERIV : Number of derivative to calculate (0<= N <=3) */
6625 /* 0 -> Position simple on base functions */
6626 /* N -> Position on base functions and derivative */
6627 /* of order 1 to N */
6629 /* OUTPUT ARGUMENTS : */
6630 /* --------------------- */
6631 /* VALBAS (NCOEFF, 0:NDERIV) : calculated value */
6633 /* d vj(t) = VALBAS(J, I) */
6637 /* IERCOD : Error code */
6639 /* 1 : Incoherence of input arguments */
6641 /* COMMONS USED : */
6642 /* -------------- */
6645 /* REFERENCES CALLED : */
6646 /* ------------------- */
6649 /* DESCRIPTION/NOTES/LIMITATIONS : */
6650 /* ----------------------------------- */
6653 /* ***********************************************************************
6656 /* ***********************************************************************
6661 /* Parameter adjustments */
6662 valbas_dim1 = *ncoeff;
6667 /* ***********************************************************************
6669 /* INITIALIZATIONS */
6670 /* ***********************************************************************
6675 /* ***********************************************************************
6678 /* ***********************************************************************
6693 iorjac = (*iordre + 1) << 1;
6695 /* (1) Generic Calculations .... */
6697 /* (1.a) Calculation of hermit polynoms */
6700 mmherm1_(moin11, &c__2, iord, hermit, &ier);
6706 /* (1.b) Evaluation of hermit polynoms */
6709 iunit = *nderiv + 1;
6710 khe = (*iordre + 1) * iunit;
6715 for (ii = 0; ii <= i__1; ++ii) {
6716 mmdrvcb_(nderiv, &c__1, &iorjac, &hermit[(ii + 3) * 6 - 18],
6717 tparam, &herm[jj - 1], &ier);
6722 mmdrvcb_(nderiv, &c__1, &iorjac, &hermit[(ii + 6) * 6 - 18],
6723 tparam, &herm[jj + khe - 1], &ier);
6733 for (ii = 0; ii <= i__1; ++ii) {
6734 AdvApp2Var_MathBase::mmpocrb_(&c__1, &iorjac, &hermit[(ii + 3) * 6 - 18], &c__1,
6735 tparam, &herm[jj - 1]);
6737 AdvApp2Var_MathBase::mmpocrb_(&c__1, &iorjac, &hermit[(ii + 6) * 6 - 18], &c__1,
6738 tparam, &herm[jj + khe - 1]);
6743 /* (1.c) Evaluation of Jacobi polynoms */
6745 ii = *ncoeff - iorjac;
6747 mmpojac_(tparam, &iorjac, &ii, nderiv, vjac, &ier);
6752 /* (1.d) Evaluation of W(t) */
6756 nwcof = advapp_max(i__1,1);
6757 AdvApp2Var_SysBase::mvriraz_(&nwcof,
6764 } else if (*iordre == 1) {
6767 } else if (*iordre == 0) {
6771 mmdrvcb_(nderiv, &c__1, &nwcof, wpoly, tparam, wval, &ier);
6776 kk1 = *ncoeff - iorjac;
6780 /* (2) Evaluation of order 0 */
6784 for (ii = 1; ii <= i__1; ++ii) {
6785 valbas[ii] = herm[jj - 1];
6790 for (ii = 1; ii <= i__1; ++ii) {
6791 valbas[ii + iorjac] = wval[0] * vjac[ii - 1];
6794 /* (3) Evaluation of order 1 */
6799 for (ii = 1; ii <= i__1; ++ii) {
6800 valbas[ii + valbas_dim1] = herm[jj - 1];
6806 for (ii = 1; ii <= i__1; ++ii) {
6807 valbas[ii + iorjac + valbas_dim1] = wval[0] * vjac[ii + kk1 - 1]
6808 + wval[1] * vjac[ii - 1];
6812 /* (4) Evaluation of order 2 */
6817 for (ii = 1; ii <= i__1; ++ii) {
6818 valbas[ii + (valbas_dim1 << 1)] = herm[jj - 1];
6823 for (ii = 1; ii <= i__1; ++ii) {
6824 valbas[ii + iorjac + (valbas_dim1 << 1)] = wval[0] * vjac[ii +
6825 kk2 - 1] + wval[1] * 2 * vjac[ii + kk1 - 1] + wval[2] *
6830 /* (5) Evaluation of order 3 */
6835 for (ii = 1; ii <= i__1; ++ii) {
6836 valbas[ii + valbas_dim1 * 3] = herm[jj - 1];
6841 for (ii = 1; ii <= i__1; ++ii) {
6842 valbas[ii + iorjac + valbas_dim1 * 3] = wval[0] * vjac[ii + kk3 -
6843 1] + wval[1] * 3 * vjac[ii + kk2 - 1] + wval[2] * 3 *
6844 vjac[ii + kk1 - 1] + wval[3] * vjac[ii - 1];
6850 /* ***********************************************************************
6852 /* ERROR PROCESSING */
6853 /* ***********************************************************************
6863 /* ***********************************************************************
6865 /* RETURN CALLING PROGRAM */
6866 /* ***********************************************************************
6872 AdvApp2Var_SysBase::maermsg_("MMPOBAS", iercod, 7L);
6877 //=======================================================================
6878 //function : AdvApp2Var_MathBase::mmpocrb_
6880 //=======================================================================
6881 int AdvApp2Var_MathBase::mmpocrb_(integer *ndimax,
6889 /* System generated locals */
6890 integer courbe_dim1, courbe_offset, i__1, i__2;
6892 /* Local variables */
6894 integer isize, nd, kcf, ncf;
6897 /* ***********************************************************************
6902 /* CALCULATE THE COORDINATES OF A POINT OF A CURVE OF GIVEN PARAMETER */
6903 /* TPARAM ( IN 2D, 3D OR MORE) */
6907 /* TOUS , MATH_ACCES :: COURBE&,PARAMETRE& , POSITIONNEMENT , &POINT
6910 /* INPUT ARGUMENTS : */
6911 /* ------------------ */
6912 /* NDIMAX : format / dimension of the curve */
6913 /* NCOEFF : Nb of coefficients of the curve */
6914 /* COURBE : Matrix of coefficients of the curve */
6915 /* NDIM : Dimension useful of the workspace */
6916 /* TPARAM : Value of the parameter where the point is calculated */
6918 /* OUTPUT ARGUMENTS : */
6919 /* ------------------- */
6920 /* PNTCRB : Coordinates of the calculated point */
6922 /* COMMONS USED : */
6923 /* ---------------- */
6927 /* REFERENCES CALLED : */
6928 /* ---------------------- */
6930 /* MIRAZ MVPSCR2 MVPSCR3 */
6932 /* DESCRIPTION/NOTES/LIMITATIONS : */
6933 /* ----------------------------------- */
6936 /* ***********************************************************************
6940 /* ***********************************************************************
6943 /* Parameter adjustments */
6944 courbe_dim1 = *ndimax;
6945 courbe_offset = courbe_dim1 + 1;
6946 courbe -= courbe_offset;
6951 AdvApp2Var_SysBase::miraz_(&isize,
6958 /* optimal processing 3d */
6960 if (*ndim == 3 && *ndimax == 3) {
6961 mvpscr3_(ncoeff, &courbe[courbe_offset], tparam, &pntcrb[1]);
6963 /* optimal processing 2d */
6965 } else if (*ndim == 2 && *ndimax == 2) {
6966 mvpscr2_(ncoeff, &courbe[courbe_offset], tparam, &pntcrb[1]);
6968 /* Any dimension - scheme of HORNER */
6970 } else if (*tparam == 0.) {
6972 for (nd = 1; nd <= i__1; ++nd) {
6973 pntcrb[nd] = courbe[nd + courbe_dim1];
6976 } else if (*tparam == 1.) {
6978 for (ncf = 1; ncf <= i__1; ++ncf) {
6980 for (nd = 1; nd <= i__2; ++nd) {
6981 pntcrb[nd] += courbe[nd + ncf * courbe_dim1];
6987 ncof2 = *ncoeff + 2;
6989 for (nd = 1; nd <= i__1; ++nd) {
6991 for (ncf = 2; ncf <= i__2; ++ncf) {
6993 pntcrb[nd] = (pntcrb[nd] + courbe[nd + kcf * courbe_dim1]) * *
6997 pntcrb[nd] += courbe[nd + courbe_dim1];
7006 //=======================================================================
7007 //function : AdvApp2Var_MathBase::mmmpocur_
7009 //=======================================================================
7010 int AdvApp2Var_MathBase::mmmpocur_(integer *ncofmx,
7018 /* System generated locals */
7019 integer courbe_dim1, courbe_offset, i__1;
7021 /* Local variables */
7026 /* ***********************************************************************
7031 /* Position of a point on curve (ncofmx,ndim). */
7035 /* TOUS , AB_SPECIFI :: COURBE&,POLYNOME&,POSITIONNEMENT,&POINT */
7037 /* INPUT ARGUMENTS : */
7038 /* ------------------ */
7039 /* NCOFMX: Format / degree of the CURVE. */
7040 /* NDIM : Dimension of the space. */
7041 /* NDEG : Degree of the polynom. */
7042 /* COURBE: Coefficients of the curve. */
7043 /* TPARAM: Parameter on the curve */
7045 /* OUTPUT ARGUMENTS : */
7046 /* ------------------- */
7047 /* TABVAL(NDIM): The resulting point (or table of values) */
7049 /* COMMONS USED : */
7050 /* ---------------- */
7052 /* REFERENCES CALLED : */
7053 /* ----------------------- */
7055 /* DESCRIPTION/NOTES/LIMITATIONS : */
7056 /* ----------------------------------- */
7059 /* ***********************************************************************
7062 /* Parameter adjustments */
7064 courbe_dim1 = *ncofmx;
7065 courbe_offset = courbe_dim1 + 1;
7066 courbe -= courbe_offset;
7071 for (nd = 1; nd <= i__1; ++nd) {
7077 for (nd = 1; nd <= i__1; ++nd) {
7078 fu = courbe[*ndeg + nd * courbe_dim1];
7079 for (i__ = *ndeg - 1; i__ >= 1; --i__) {
7080 fu = fu * *tparam + courbe[i__ + nd * courbe_dim1];
7090 //=======================================================================
7091 //function : mmpojac_
7093 //=======================================================================
7094 int mmpojac_(doublereal *tparam,
7104 /* Initialized data */
7108 /* System generated locals */
7109 integer valjac_dim1, i__1, i__2;
7111 /* Local variables */
7112 doublereal cofa, cofb, denom, tnorm[100];
7113 integer ii, jj, kk1, kk2;
7114 doublereal aux1, aux2;
7117 /* ***********************************************************************
7122 /* Positioning on Jacobi polynoms and their derivatives */
7123 /* successive by a recurrent algorithm */
7127 /* RESERVE, POSITIONING, JACOBI */
7129 /* INPUT ARGUMENTS : */
7130 /* -------------------- */
7131 /* TPARAM : Parameter for which positioning is done. */
7132 /* IORDRE : Order of hermit-?? (-1,0,1, or 2) */
7133 /* NCOEFF : Number of coeeficients of polynoms (Nb of value to */
7135 /* NDERIV : Number of derivative to calculate (0<= N <=3) */
7136 /* 0 -> Position simple on jacobi functions */
7137 /* N -> Position on jacobi functions and their */
7138 /* derivatives of order 1 to N. */
7140 /* OUTPUT ARGUMENTS : */
7141 /* --------------------- */
7142 /* VALJAC (NCOEFF, 0:NDERIV) : the calculated values */
7144 /* d vj(t) = VALJAC(J, I) */
7148 /* IERCOD : Error Code */
7150 /* 1 : Incoherence of input arguments */
7152 /* COMMONS USED : */
7153 /* ------------------ */
7156 /* REFERENCES CALLED : */
7157 /* --------------------- */
7160 /* DESCRIPTION/NOTES/LIMITATIONS : */
7161 /* ----------------------------------- */
7164 /* ***********************************************************************
7167 /* ***********************************************************************
7171 /* static varaibles */
7175 /* Parameter adjustments */
7176 valjac_dim1 = *ncoeff;
7181 /* ***********************************************************************
7183 /* INITIALISATIONS */
7184 /* ***********************************************************************
7189 /* ***********************************************************************
7192 /* ***********************************************************************
7198 if (*ncoeff > 100) {
7202 /* --- Calculation of norms */
7204 /* IF (NCOEFF.GT.NBCOF) THEN */
7206 for (ii = 1; ii <= i__1; ++ii) {
7210 for (jj = 1; jj <= i__2; ++jj) {
7211 aux2 = aux2 * (doublereal) (kk1 + *iordre + jj) / (doublereal) (
7214 i__2 = (*iordre << 1) + 1;
7215 tnorm[ii - 1] = sqrt(aux2 * (kk1 * 2. + (*iordre << 1) + 1) / pow__ii(&
7223 /* --- Trivial Positions ----- */
7226 aux1 = (doublereal) (*iordre + 1);
7227 valjac[2] = aux1 * *tparam;
7230 valjac[valjac_dim1 + 1] = 0.;
7231 valjac[valjac_dim1 + 2] = aux1;
7234 valjac[(valjac_dim1 << 1) + 1] = 0.;
7235 valjac[(valjac_dim1 << 1) + 2] = 0.;
7238 valjac[valjac_dim1 * 3 + 1] = 0.;
7239 valjac[valjac_dim1 * 3 + 2] = 0.;
7244 /* --- Positioning by recurrence */
7247 for (ii = 3; ii <= i__1; ++ii) {
7251 aux1 = (doublereal) (*iordre + kk2);
7253 cofa = aux2 * (aux2 + 1) * (aux2 + 2);
7254 cofb = (aux2 + 2) * -2. * aux1 * aux1;
7255 denom = kk1 * 2. * (kk2 + (*iordre << 1) + 1) * aux2;
7259 valjac[ii] = (cofa * *tparam * valjac[kk1] + cofb * valjac[kk2]) *
7263 valjac[ii + valjac_dim1] = (cofa * *tparam * valjac[kk1 +
7264 valjac_dim1] + cofa * valjac[kk1] + cofb * valjac[kk2 +
7265 valjac_dim1]) * denom;
7268 valjac[ii + (valjac_dim1 << 1)] = (cofa * *tparam * valjac[
7269 kk1 + (valjac_dim1 << 1)] + cofa * 2 * valjac[kk1 +
7270 valjac_dim1] + cofb * valjac[kk2 + (valjac_dim1 << 1)]
7275 valjac[ii + valjac_dim1 * 3] = (cofa * *tparam * valjac[kk1 +
7276 valjac_dim1 * 3] + cofa * 3 * valjac[kk1 + (
7277 valjac_dim1 << 1)] + cofb * valjac[kk2 + valjac_dim1 *
7283 /* ---> Normalization */
7286 for (ii = 1; ii <= i__1; ++ii) {
7288 for (jj = 0; jj <= i__2; ++jj) {
7289 valjac[ii + jj * valjac_dim1] = tnorm[ii - 1] * valjac[ii + jj *
7296 /* ***********************************************************************
7298 /* PROCESSING OF ERRORS */
7299 /* ***********************************************************************
7307 /* ***********************************************************************
7309 /* RETURN CALLING PROGRAM */
7310 /* ***********************************************************************
7316 AdvApp2Var_SysBase::maermsg_("MMPOJAC", iercod, 7L);
7321 //=======================================================================
7322 //function : AdvApp2Var_MathBase::mmposui_
7324 //=======================================================================
7325 int AdvApp2Var_MathBase::mmposui_(integer *dimmat,
7332 /* System generated locals */
7335 /* Local variables */
7337 integer imin, jmin, i__, j, k;
7340 /* ***********************************************************************
7345 /* FILL THE TABLE OF POSITIONING POSUIV WHICH ALLOWS TO */
7346 /* PARSE BY COLUMN THE INFERIOR TRIANGULAR PART OF THE */
7347 /* MATRIX IN FORM OF PROFILE */
7352 /* RESERVE, MATRIX, PROFILE */
7354 /* INPUT ARGUMENTS : */
7355 /* -------------------- */
7357 /* NISTOC: NUMBER OF COEFFICIENTS IN THE PROFILE */
7358 /* DIMMAT: NUMBER OF LINE OF THE SYMMETRIC SQUARE MATRIX */
7359 /* APOSIT: TABLE OF POSITIONING OF STORAGE TERMS */
7360 /* APOSIT(1,I) CONTAINS THE NUMBER OF TERMES-1 ON LINE
7361 /* I IN THE PROFILE OF THE MATRIX */
7362 /* APOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF DIAGONAL TERM
7366 /* OUTPUT ARGUMENTS : */
7367 /* --------------------- */
7368 /* POSUIV: POSUIV(K) (WHERE K IS THE INDEX OF STORAGE OF MAT(I,J)) */
7369 /* CONTAINS THE SMALLEST NUMBER IMIN>I OF THE LINE THAT */
7370 /* POSSESSES A TERM MAT(IMIN,J) THAT IS IN THE PROFILE. */
7371 /* IF THERE IS NO TERM MAT(IMIN,J) IN THE PROFILE THEN POSUIV(K)=-1 */
7374 /* COMMONS USED : */
7375 /* ------------------ */
7378 /* REFERENCES CALLED : */
7379 /* --------------------- */
7382 /* DESCRIPTION/NOTES/LIMITATIONS : */
7383 /* ----------------------------------- */
7386 /* ***********************************************************************
7389 /* ***********************************************************************
7394 /* ***********************************************************************
7396 /* INITIALIZATIONS */
7397 /* ***********************************************************************
7400 /* Parameter adjustments */
7405 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
7407 AdvApp2Var_SysBase::mgenmsg_("MMPOSUI", 7L);
7412 /* ***********************************************************************
7415 /* ***********************************************************************
7421 for (i__ = 1; i__ <= i__1; ++i__) {
7422 jmin = i__ - aposit[(i__ << 1) + 1];
7424 for (j = jmin; j <= i__2; ++j) {
7427 while(! trouve && imin <= *dimmat) {
7428 if (imin - aposit[(imin << 1) + 1] <= j) {
7434 k = aposit[(i__ << 1) + 2] - i__ + j;
7449 /* ***********************************************************************
7451 /* ERROR PROCESSING */
7452 /* ***********************************************************************
7458 /* ***********************************************************************
7460 /* RETURN CALLING PROGRAM */
7461 /* ***********************************************************************
7466 /* ___ DESALLOCATION, ... */
7468 AdvApp2Var_SysBase::maermsg_("MMPOSUI", iercod, 7L);
7470 AdvApp2Var_SysBase::mgsomsg_("MMPOSUI", 7L);
7475 //=======================================================================
7476 //function : AdvApp2Var_MathBase::mmresol_
7478 //=======================================================================
7479 int AdvApp2Var_MathBase::mmresol_(integer *hdimen,
7497 integer c__100 = 100;
7499 /* System generated locals */
7502 /* Local variables */
7504 doublereal* mcho = 0;
7505 integer jmin, jmax, i__, j, k, l;
7506 intptr_t iofv1, iofv2, iofv3, iofv4;
7507 doublereal *v1 = 0, *v2 = 0, *v3 = 0, *v4 = 0;
7508 integer deblig, dimhch;
7509 doublereal* hchole = 0;
7510 intptr_t iofmch, iofmam, iofhch;
7511 doublereal* matsym = 0;
7517 /* ***********************************************************************
7522 /* SOLUTION OF THE SYSTEM */
7529 /* RESERVE, SOLUTION, SYSTEM, LAGRANGIAN */
7531 /* INPUT ARGUMENTS : */
7532 /* -------------------- */
7533 /* HDIMEN: NOMBER OF LINE (OR COLUMN) OF THE HESSIAN MATRIX */
7534 /* GDIMEN: NOMBER OF LINE OF THE MATRIX OF CONSTRAINTS */
7535 /* HNSTOC: NOMBErS OF TERMS IN THE PROFILE OF HESSIAN MATRIX
7537 /* GNSTOC: NOMBERS OF TERMS IN THE PROFILE OF THE MATRIX OF CONSTRAINTS */
7538 /* MNSTOC: NOMBERS OF TERMS IN THE PROFILE OF THE MATRIX M= G H t(G) */
7539 /* where H IS THE HESSIAN MATRIX AND G IS THE MATRIX OF CONSTRAINTS */
7540 /* MATSYH: TRIANGULAR INFERIOR PART OF THE HESSIAN MATRIX
7541 /* IN FORM OF PROFILE */
7542 /* MATSYG: MATRIX OF CONSTRAINTS IN FORM OF PROFILE */
7543 /* VECSYH: VECTOR OF THE SECOND MEMBER ASSOCIATED TO MATSYH */
7544 /* VECSYG: VECTOR OF THE SECOND MEMBER ASSOCIATED TO MATSYG */
7545 /* HPOSIT: TABLE OF POSITIONING OF THE HESSIAN MATRIX */
7546 /* HPOSIT(1,I) CONTAINS THE NUMBER OF TERMS -1 */
7547 /* WHICH ARE IN THE PROFILE AT LINE I */
7548 /* HPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF TERM */
7549 /* DIAGONAL OF THE MATRIX AT LINE I */
7550 /* HPOSUI: TABLE ALLOWING TO PARSE THE HESSIAN MATRIX BY COLUMN */
7551 /* IN FORM OF PROFILE */
7552 /* HPOSUI(K) CONTAINS THE NUMBER OF LINE IMIN FOLLOWING THE CURRENT LINE*/
7553 /* I WHERE H(I,J)=MATSYH(K) AS IT EXISTS IN THE */
7554 /* SAME COLUMN J A TERM IN THE PROFILE OF LINE IMIN */
7555 /* IF SUCH TERM DOES NOT EXIST IMIN=-1 */
7556 /* GPOSIT: TABLE OF POSITIONING OF THE MATRIX OF CONSTRAINTS */
7557 /* GPOSIT(1,I) CONTAINS THE NUMBER OF TERMS OF LINE I */
7558 /* WHICH ARE IN THE PROFILE */
7559 /* GPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF THE LAST TERM
7560 /* OF LINE I WHICH IS IN THE PROFILE */
7561 /* GPOSIT(3,I) CONTAINS THE NUMBER OF COLUMN CORRESPONDING */
7562 /* TO THE FIRST TERM OF LINE I WHICH IS IN THE PROFILE */
7563 /* MMPOSUI, MPOSIT: SAME STRUCTURE AS HPOSUI, BUT FOR MATRIX
7567 /* OUTPUT ARGUMENTS : */
7568 /* --------------------- */
7569 /* VECSOL: VECTOR SOLUTION V OF THE SYSTEM */
7570 /* IERCOD: ERROR CODE */
7572 /* COMMONS USED : */
7573 /* ------------------ */
7576 /* REFERENCES CALLED : */
7577 /* --------------------- */
7580 /* DESCRIPTION/NOTES/LIMITATIONS : */
7581 /* ----------------------------------- */
7583 /* ***********************************************************************
7586 /* ***********************************************************************
7589 /* ***********************************************************************
7591 /* INITIALISATIONS */
7592 /* ***********************************************************************
7595 /* Parameter adjustments */
7608 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
7610 AdvApp2Var_SysBase::mgenmsg_("MMRESOL", 7L);
7621 /* ***********************************************************************
7624 /* ***********************************************************************
7627 /* Dynamic allocation */
7628 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
7629 anAdvApp2Var_SysBase.macrar8_(hdimen, &c__100, v1, &iofv1, &ier);
7633 dimhch = hposit[(*hdimen << 1) + 2];
7634 anAdvApp2Var_SysBase.macrar8_(&dimhch, &c__100, hchole, &iofhch, &ier);
7639 /* solution of system 1 H V1 = b */
7640 /* where H=MATSYH and b=VECSYH */
7642 mmchole_(hnstoc, hdimen, &matsyh[1], &hposit[3], &hposui[1], &hchole[
7647 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1], &vecsyh[
7648 1], &v1[iofv1], &ier);
7653 /* Case when there are constraints */
7657 /* Calculate the vector of the second member V2=G H(-1) b -c = G v1-c */
7658 /* of system of unknown Lagrangian vector MULTIP */
7659 /* where G=MATSYG */
7662 anAdvApp2Var_SysBase.macrar8_(gdimen, &c__100, v2, &iofv2, &ier);
7666 anAdvApp2Var_SysBase.macrar8_(hdimen, &c__100, v3, &iofv3, &ier);
7670 anAdvApp2Var_SysBase.macrar8_(gdimen, &c__100, v4, &iofv4, &ier);
7674 anAdvApp2Var_SysBase.macrar8_(mnstoc, &c__100, matsym, &iofmam, &ier);
7680 mmatvec_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v1[iofv1], &
7681 deblig, &v2[iofv2], &ier);
7686 for (i__ = 1; i__ <= i__1; ++i__) {
7687 v2[i__ + iofv2 - 1] -= vecsyg[i__];
7690 /* Calculate the matrix M= G H(-1) t(G) */
7691 /* RESOL DU SYST 2 : H qi = gi */
7692 /* where is a vector column of t(G) */
7694 /* then calculate G qi */
7695 /* then construct M in form of profile */
7700 for (i__ = 1; i__ <= i__1; ++i__) {
7701 AdvApp2Var_SysBase::mvriraz_(hdimen, &v1[iofv1]);
7702 AdvApp2Var_SysBase::mvriraz_(hdimen, &v3[iofv3]);
7703 AdvApp2Var_SysBase::mvriraz_(gdimen, &v4[iofv4]);
7704 jmin = gposit[i__ * 3 + 3];
7705 jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
7706 aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
7708 for (j = jmin; j <= i__2; ++j) {
7710 v1[j + iofv1 - 1] = matsyg[k];
7712 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1],
7713 &v1[iofv1], &v3[iofv3], &ier);
7719 mmatvec_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v3[
7720 iofv3], &deblig, &v4[iofv4], &ier);
7725 k = mposit[(i__ << 1) + 2];
7726 matsym[k + iofmam - 1] = v4[i__ + iofv4 - 1];
7727 while(mmposui[k] > 0) {
7729 k = mposit[(l << 1) + 2] - l + i__;
7730 matsym[k + iofmam - 1] = v4[l + iofv4 - 1];
7735 /* SOLVE SYST 3 M L = V2 */
7739 AdvApp2Var_SysBase::mvriraz_(gdimen, &v4[iofv4]);
7740 anAdvApp2Var_SysBase.macrar8_(mnstoc, &c__100, mcho, &iofmch, &ier);
7744 mmchole_(mnstoc, gdimen, &matsym[iofmam], &mposit[3], &mmposui[1], &
7745 mcho[iofmch], &ier);
7749 mmrslss_(mnstoc, gdimen, &mcho[iofmch], &mposit[3], &mmposui[1], &v2[
7750 iofv2], &v4[iofv4], &ier);
7756 /* CALCULATE THE VECTOR OF THE SECOND MEMBER OF THE SYSTEM Hx = b - t(G) L
7760 AdvApp2Var_SysBase::mvriraz_(hdimen, &v1[iofv1]);
7761 mmtmave_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v4[iofv4], &
7767 for (i__ = 1; i__ <= i__1; ++i__) {
7768 v1[i__ + iofv1 - 1] = vecsyh[i__] - v1[i__ + iofv1 - 1];
7771 /* RESOL SYST 4 Hx = b - t(G) L */
7774 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1], &v1[
7775 iofv1], &vecsol[1], &ier);
7781 for (i__ = 1; i__ <= i__1; ++i__) {
7782 vecsol[i__] = v1[i__ + iofv1 - 1];
7788 /* ***********************************************************************
7790 /* PROCESSING OF ERRORS */
7791 /* ***********************************************************************
7800 AdvApp2Var_SysBase::mswrdbg_("MMRESOL : PROBLEM WITH DIMMAT", 30L);
7803 /* ***********************************************************************
7805 /* RETURN CALLING PROGRAM */
7806 /* ***********************************************************************
7811 /* ___ DESALLOCATION, ... */
7812 anAdvApp2Var_SysBase.macrdr8_(hdimen, &c__100, v1, &iofv1, &ier);
7813 if (*iercod == 0 && ier > 0) {
7816 anAdvApp2Var_SysBase.macrdr8_(&dimhch, &c__100, hchole, &iofhch, &ier);
7817 if (*iercod == 0 && ier > 0) {
7820 anAdvApp2Var_SysBase.macrdr8_(gdimen, &c__100, v2, &iofv2, &ier);
7821 if (*iercod == 0 && ier > 0) {
7824 anAdvApp2Var_SysBase.macrdr8_(hdimen, &c__100, v3, &iofv3, &ier);
7825 if (*iercod == 0 && ier > 0) {
7828 anAdvApp2Var_SysBase.macrdr8_(gdimen, &c__100, v4, &iofv4, &ier);
7829 if (*iercod == 0 && ier > 0) {
7832 anAdvApp2Var_SysBase.macrdr8_(mnstoc, &c__100, matsym, &iofmam, &ier);
7833 if (*iercod == 0 && ier > 0) {
7836 anAdvApp2Var_SysBase.macrdr8_(mnstoc, &c__100, mcho, &iofmch, &ier);
7837 if (*iercod == 0 && ier > 0) {
7841 AdvApp2Var_SysBase::maermsg_("MMRESOL", iercod, 7L);
7843 AdvApp2Var_SysBase::mgsomsg_("MMRESOL", 7L);
7848 //=======================================================================
7849 //function : mmrslss_
7851 //=======================================================================
7852 int mmrslss_(integer *,//mxcoef,
7857 doublereal *mscnmbr,
7861 /* System generated locals */
7864 /* Local variables */
7868 integer pointe, ptcour;
7870 /* ***********************************************************************
7875 /* Solves linear system SS x = b where S is a */
7876 /* triangular lower matrix given in form of profile */
7880 /* RESERVE, MATRICE_PROFILE, RESOLUTION, CHOLESKI */
7882 /* INPUT ARGUMENTS : */
7883 /* -------------------- */
7884 /* MXCOEF : Maximum number of non-null coefficient in the matrix */
7885 /* DIMENS : Dimension of the matrix */
7886 /* SMATRI(MXCOEF) : Values of coefficients of the matrix */
7887 /* SPOSIT(2,DIMENS): */
7888 /* SPOSIT(1,*) : Distance diagonal-extremity of the line */
7889 /* SPOSIT(2,*) : Position of diagonal terms in AMATRI */
7890 /* POSUIV(MXCOEF): first line inferior not out of profile */
7891 /* MSCNMBR(DIMENS): Vector second member of the equation */
7893 /* OUTPUT ARGUMENTS : */
7894 /* --------------------- */
7895 /* SOLUTI(NDIMEN) : Result vector */
7896 /* IERCOD : Error code 0 : ok */
7898 /* COMMONS USED : */
7899 /* ------------------ */
7902 /* REFERENCES CALLED : */
7903 /* --------------------- */
7906 /* DESCRIPTION/NOTES/LIMITATIONS : */
7907 /* ----------------------------------- */
7909 /* SS is the decomposition of choleski of a symmetric matrix */
7910 /* defined postive, that can result from routine MMCHOLE. */
7912 /* For a full matrix it is possible to use MRSLMSC */
7914 /* LEVEL OF DEBUG = 4 */
7916 /* ***********************************************************************
7919 /* ***********************************************************************
7924 /* ***********************************************************************
7926 /* INITIALISATIONS */
7927 /* ***********************************************************************
7930 /* Parameter adjustments */
7938 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
7940 AdvApp2Var_SysBase::mgenmsg_("MMRSLSS", 7L);
7944 /* ***********************************************************************
7947 /* ***********************************************************************
7950 /* ----- Solution of Sw = b */
7953 for (i__ = 1; i__ <= i__1; ++i__) {
7955 pointe = sposit[(i__ << 1) + 2];
7958 for (j = i__ - sposit[(i__ << 1) + 1]; j <= i__2; ++j) {
7959 somme += smatri[pointe - (i__ - j)] * soluti[j];
7962 soluti[i__] = (mscnmbr[i__] - somme) / smatri[pointe];
7965 /* ----- Solution of S u = w */
7967 for (i__ = *dimens; i__ >= 1; --i__) {
7969 pointe = sposit[(i__ << 1) + 2];
7973 ptcour = sposit[(j << 1) + 2] - (j - i__);
7974 somme += smatri[ptcour] * soluti[j];
7978 soluti[i__] = (soluti[i__] - somme) / smatri[pointe];
7983 /* ***********************************************************************
7985 /* ERROR PROCESSING */
7986 /* ***********************************************************************
7990 /* ***********************************************************************
7992 /* RETURN PROGRAM CALLING */
7993 /* ***********************************************************************
7998 AdvApp2Var_SysBase::maermsg_("MMRSLSS", iercod, 7L);
8000 AdvApp2Var_SysBase::mgsomsg_("MMRSLSS", 7L);
8005 //=======================================================================
8006 //function : mmrslw_
8008 //=======================================================================
8009 int mmrslw_(integer *normax,
8017 /* System generated locals */
8018 integer abmatr_dim1, abmatr_offset, xmatri_dim1, xmatri_offset, i__1,
8022 /* Local variables */
8029 /* **********************************************************************
8034 /* Solution of a linear system A.x = B of N equations to N */
8035 /* unknown by Gauss method (partial pivot) or : */
8036 /* A is matrix NORDRE * NORDRE, */
8037 /* B is matrix NORDRE (lines) * NDIMEN (columns), */
8038 /* x is matrix NORDRE (lines) * NDIMEN (columns). */
8039 /* In this program, A and B are stored in matrix ABMATR */
8040 /* the lines and columns which of were inverted. ABMATR(k,j) is */
8041 /* term A(j,k) if k <= NORDRE, B(j,k-NORDRE) otherwise (see example). */
8045 /* TOUS, MATH_ACCES::EQUATION&, MATRICE&, RESOLUTION, GAUSS, &SOLUTION */
8047 /* INPUT ARGUMENTS : */
8048 /* ------------------ */
8049 /* NORMAX : Max size of the first index of XMATRI. This argument */
8050 /* serves only for the declaration of dimension of XMATRI and should be */
8051 /* above or equal to NORDRE. */
8052 /* NORDRE : Order of the matrix i.e. number of equations and */
8053 /* unknown quantities of the linear system to be solved. */
8054 /* NDIMEN : Number of the second member. */
8055 /* EPSPIV : Minimal value of a pivot. If during the calculation */
8056 /* the absolute value of the pivot is below EPSPIV, the */
8057 /* system of equations is declared singular. EPSPIV should */
8058 /* be a "small" real. */
8060 /* ABMATR(NORDRE+NDIMEN,NORDRE) : Auxiliary matrix containing */
8061 /* matrix A and matrix B. */
8063 /* OUTPUT ARGUMENTS : */
8064 /* ------------------- */
8065 /* XMATRI : Matrix containing NORDRE*NDIMEN solutions. */
8066 /* IERCOD=0 shows that all solutions are calculated. */
8067 /* IERCOD=1 shows that the matrix is of lower rank than NORDRE */
8068 /* (the system is singular). */
8070 /* COMMONS USED : */
8071 /* ---------------- */
8073 /* REFERENCES CALLED : */
8074 /* ----------------------- */
8076 /* DESCRIPTION/NOTES/LIMITATIONS : */
8077 /* ----------------------------------- */
8078 /* ATTENTION : the indices of line and column are inverted */
8079 /* compared to usual indices. */
8081 /* a1*x + b1*y = c1 */
8082 /* a2*x + b2*y = c2 */
8083 /* should be represented by matrix ABMATR : */
8085 /* ABMATR(1,1) = a1 ABMATR(1,2) = a2 */
8086 /* ABMATR(2,1) = b1 ABMATR(2,2) = b2 */
8087 /* ABMATR(3,1) = c1 ABMATR(3,2) = c2 */
8089 /* To solve this system, it is necessary to set : */
8091 /* NORDRE = 2 (there are 2 equations with 2 unknown values), */
8092 /* NDIMEN = 1 (there is only one second member), */
8093 /* any NORMAX can be taken >= NORDRE. */
8095 /* To use this routine, it is recommended to use one of */
8096 /* interfaces : MMRSLWI or MMMRSLWD. */
8098 /* **********************************************************************
8101 /* Name of the routine */
8103 /* INTEGER IBB,MNFNDEB */
8106 /* IF (IBB.GE.2) CALL MGENMSG(NOMPR) */
8107 /* Parameter adjustments */
8108 xmatri_dim1 = *normax;
8109 xmatri_offset = xmatri_dim1 + 1;
8110 xmatri -= xmatri_offset;
8111 abmatr_dim1 = *nordre + *ndimen;
8112 abmatr_offset = abmatr_dim1 + 1;
8113 abmatr -= abmatr_offset;
8118 /* *********************************************************************
8120 /* Triangulation of matrix ABMATR. */
8121 /* *********************************************************************
8125 for (kk = 1; kk <= i__1; ++kk) {
8127 /* ---------- Find max pivot in column KK. ------------
8133 for (jj = kk; jj <= i__2; ++jj) {
8134 akj = (d__1 = abmatr[kk + jj * abmatr_dim1], advapp_abs(d__1));
8145 /* --------- Swapping of line KPIV with line KK. ------
8149 i__2 = *nordre + *ndimen;
8150 for (jj = kk; jj <= i__2; ++jj) {
8151 akj = abmatr[jj + kk * abmatr_dim1];
8152 abmatr[jj + kk * abmatr_dim1] = abmatr[jj + kpiv *
8154 abmatr[jj + kpiv * abmatr_dim1] = akj;
8159 /* ---------- Removal and triangularization. -----------
8162 pivot = -abmatr[kk + kk * abmatr_dim1];
8164 for (ii = kk + 1; ii <= i__2; ++ii) {
8165 akj = abmatr[kk + ii * abmatr_dim1] / pivot;
8166 i__3 = *nordre + *ndimen;
8167 for (jj = kk + 1; jj <= i__3; ++jj) {
8168 abmatr[jj + ii * abmatr_dim1] += akj * abmatr[jj + kk *
8179 /* *********************************************************************
8181 /* Solution of the system of triangular equations. */
8182 /* Matrix ABMATR(NORDRE+JJ,II), contains second members */
8183 /* of the system for 1<=j<=NDIMEN and 1<=i<=NORDRE. */
8184 /* *********************************************************************
8188 /* ---------------- Calculation of solutions by ascending. -----------------
8191 for (kk = *nordre; kk >= 1; --kk) {
8192 pivot = abmatr[kk + kk * abmatr_dim1];
8194 for (ii = 1; ii <= i__1; ++ii) {
8195 akj = abmatr[ii + *nordre + kk * abmatr_dim1];
8197 for (jj = kk + 1; jj <= i__2; ++jj) {
8198 akj -= abmatr[jj + kk * abmatr_dim1] * xmatri[jj + ii *
8202 xmatri[kk + ii * xmatri_dim1] = akj / pivot;
8209 /* ------If the absolute value of a pivot is smaller than --------
8210 /* ---------- EPSPIV: return the code of error. ------------
8220 AdvApp2Var_SysBase::maermsg_("MMRSLW ", iercod, 7L);
8222 /* IF (IBB.GE.2) CALL MGSOMSG(NOMPR) */
8226 //=======================================================================
8227 //function : AdvApp2Var_MathBase::mmmrslwd_
8229 //=======================================================================
8230 int AdvApp2Var_MathBase::mmmrslwd_(integer *normax,
8241 /* System generated locals */
8242 integer amat_dim1, amat_offset, bmat_dim1, bmat_offset, xmat_dim1,
8243 xmat_offset, aaux_dim1, aaux_offset, i__1, i__2;
8245 /* Local variables */
8249 /* IMPLICIT DOUBLE PRECISION (A-H,O-Z) */
8250 /* IMPLICIT INTEGER (I-N) */
8253 /* **********************************************************************
8258 /* Solution of a linear system by Gauss method where */
8259 /* the second member is a table of vectors. Method of partial pivot. */
8263 /* ALL, MATH_ACCES :: */
8264 /* SYSTEME&,EQUATION&, RESOLUTION,GAUSS ,&VECTEUR */
8266 /* INPUT ARGUMENTS : */
8267 /* ------------------ */
8268 /* NORMAX : Max. Dimension of AMAT. */
8269 /* NORDRE : Order of the matrix. */
8270 /* NDIM : Number of columns of BMAT and XMAT. */
8271 /* AMAT(NORMAX,NORDRE) : The processed matrix. */
8272 /* BMAT(NORMAX,NDIM) : The matrix of second member. */
8273 /* XMAT(NORMAX,NDIM) : The matrix of solutions. */
8274 /* EPSPIV : Min value of a pivot. */
8276 /* OUTPUT ARGUMENTS : */
8277 /* ------------------- */
8278 /* AAUX(NORDRE+NDIM,NORDRE) : Auxiliary matrix. */
8279 /* XMAT(NORMAX,NDIM) : Matrix of solutions. */
8280 /* IERCOD=0 shows that solutions in XMAT are valid. */
8281 /* IERCOD=1 shows that matrix AMAT is of lower rank than NORDRE. */
8283 /* COMMONS USED : */
8284 /* ---------------- */
8288 /* REFERENCES CALLED : */
8289 /* ---------------------- */
8291 /* MAERMSG MGENMSG MGSOMSG */
8292 /* MMRSLW I*4 MNFNDEB */
8294 /* DESCRIPTION/NOTES/LIMITATIONS : */
8295 /* ----------------------------------- */
8296 /* ATTENTION : lines and columns are located in usual order : */
8297 /* 1st index = index line */
8298 /* 2nd index = index column */
8299 /* Example, the system : */
8300 /* a1*x + b1*y = c1 */
8301 /* a2*x + b2*y = c2 */
8302 /* is represented by matrix AMAT : */
8304 /* AMAT(1,1) = a1 AMAT(2,1) = a2 */
8305 /* AMAT(1,2) = b1 AMAT(2,2) = b2 */
8307 /* The first index is the index of line, the second index */
8308 /* is the index of columns (Compare with MMRSLWI which is faster). */
8311 /* **********************************************************************
8314 /* Name of the routine */
8316 /* Parameter adjustments */
8317 amat_dim1 = *normax;
8318 amat_offset = amat_dim1 + 1;
8319 amat -= amat_offset;
8320 xmat_dim1 = *normax;
8321 xmat_offset = xmat_dim1 + 1;
8322 xmat -= xmat_offset;
8323 aaux_dim1 = *nordre + *ndim;
8324 aaux_offset = aaux_dim1 + 1;
8325 aaux -= aaux_offset;
8326 bmat_dim1 = *normax;
8327 bmat_offset = bmat_dim1 + 1;
8328 bmat -= bmat_offset;
8331 ibb = AdvApp2Var_SysBase::mnfndeb_();
8333 AdvApp2Var_SysBase::mgenmsg_("MMMRSLW", 7L);
8336 /* Initialization of the auxiliary matrix. */
8339 for (i__ = 1; i__ <= i__1; ++i__) {
8341 for (j = 1; j <= i__2; ++j) {
8342 aaux[j + i__ * aaux_dim1] = amat[i__ + j * amat_dim1];
8348 /* Second member. */
8351 for (i__ = 1; i__ <= i__1; ++i__) {
8353 for (j = 1; j <= i__2; ++j) {
8354 aaux[j + *nordre + i__ * aaux_dim1] = bmat[i__ + j * bmat_dim1];
8360 /* Solution of the system of equations. */
8362 mmrslw_(normax, nordre, ndim, epspiv, &aaux[aaux_offset], &xmat[
8363 xmat_offset], iercod);
8367 AdvApp2Var_SysBase::maermsg_("MMMRSLW", iercod, 7L);
8370 AdvApp2Var_SysBase::mgsomsg_("MMMRSLW", 7L);
8375 //=======================================================================
8376 //function : AdvApp2Var_MathBase::mmrtptt_
8378 //=======================================================================
8379 int AdvApp2Var_MathBase::mmrtptt_(integer *ndglgd,
8383 integer ideb, nmod2, nsur2, ilong, ibb;
8386 /* **********************************************************************
8391 /* Extracts from Common LDGRTL the STRICTLY positive roots of the */
8392 /* Legendre polynom of degree NDGLGD, for 2 <= NDGLGD <= 61. */
8396 /* TOUS, AB_SPECIFI::COMMON&, EXTRACTION, &RACINE, &LEGENDRE. */
8398 /* INPUT ARGUMENTS : */
8399 /* ------------------ */
8400 /* NDGLGD : Mathematic degree of Legendre polynom. */
8401 /* This degree should be above or equal to 2 and */
8402 /* below or equal to 61. */
8404 /* OUTPUT ARGUMENTS : */
8405 /* ------------------- */
8406 /* RTLEGD : The table of strictly positive roots of */
8407 /* Legendre polynom of degree NDGLGD. */
8409 /* COMMONS USED : */
8410 /* ---------------- */
8412 /* REFERENCES CALLED : */
8413 /* ----------------------- */
8415 /* DESCRIPTION/NOTES/LIMITATIONS : */
8416 /* ----------------------------------- */
8417 /* ATTENTION: the condition on NDEGRE ( 2 <= NDEGRE <= 61) is not */
8418 /* tested. The caller should make the test. */
8421 /* **********************************************************************
8423 /* Nome of the routine */
8426 /* Common MLGDRTL: */
8427 /* This common includes POSITIVE roots of Legendre polynoms */
8428 /* AND the weight of Gauss quadrature formulas on all */
8429 /* POSITIVE roots of Legendre polynoms. */
8432 /* ***********************************************************************
8437 /* The common of Legendre roots. */
8443 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
8444 /* ----------------------------------- */
8447 /* ***********************************************************************
8453 /* ROOTAB : Table of all rotts of Legendre polynoms */
8454 /* between [0,1]. They are ranked for degrees increasing from 2 to 61. */
8455 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
8456 /* The address is the same. */
8457 /* HI0TAB : Table of Legendre interpolators for root x=0 */
8458 /* the polynoms of UNEVEN degree. */
8459 /* RTLTB0 : Table of Li(uk) where uk are roots of a */
8460 /* Legendre polynom of EVEN degree. */
8461 /* RTLTB1 : Table of Li(uk) where uk are roots of a */
8462 /* Legendre polynom of UNEVEN degree. */
8465 /************************************************************************
8467 /* Parameter adjustments */
8471 ibb = AdvApp2Var_SysBase::mnfndeb_();
8473 AdvApp2Var_SysBase::mgenmsg_("MMRTPTT", 7L);
8479 nsur2 = *ndglgd / 2;
8480 nmod2 = *ndglgd % 2;
8483 ideb = nsur2 * (nsur2 - 1) / 2 + 1;
8484 AdvApp2Var_SysBase::mcrfill_(&ilong,
8485 &mlgdrtl_.rootab[ideb + nmod2 * 465 - 1],
8488 /* ----------------------------- The end --------------------------------
8493 AdvApp2Var_SysBase::mgsomsg_("MMRTPTT", 7L);
8498 //=======================================================================
8499 //function : AdvApp2Var_MathBase::mmsrre2_
8501 //=======================================================================
8502 int AdvApp2Var_MathBase::mmsrre2_(doublereal *tparam,
8510 /* System generated locals */
8513 /* Local variables */
8514 integer ideb, ifin, imil, ibb;
8516 /* ***********************************************************************
8522 /* Find the interval corresponding to a valueb given in */
8523 /* increasing order of real numbers with double precision. */
8527 /* TOUS,MATH_ACCES::TABLEAU&,POINT&,CORRESPONDANCE,&RANG */
8529 /* INPUT ARGUMENTS : */
8530 /* ------------------ */
8532 /* TPARAM : Value to be tested. */
8533 /* NBRVAL : Size of TABLEV */
8534 /* TABLEV : Table of reals. */
8535 /* EPSIL : Epsilon of precision */
8537 /* OUTPUT ARGUMENTS : */
8538 /* ------------------- */
8540 /* NUMINT : Number of the interval (between 1 and NBRVAL-1). */
8541 /* ITYPEN : = 0 TPARAM is inside the interval NUMINT */
8542 /* = 1 : TPARAM corresponds to the lower limit of */
8543 /* the provided interval. */
8544 /* = 2 : TPARAM corresponds to the upper limit of */
8545 /* the provided interval. */
8547 /* IERCOD : Error code. */
8549 /* = 1 : TABLEV does not contain enough elements. */
8550 /* = 2 : TPARAM out of limits of TABLEV. */
8552 /* COMMONS USED : */
8553 /* ---------------- */
8555 /* REFERENCES CALLED : */
8556 /* ------------------- */
8558 /* DESCRIPTION/NOTES/LIMITATIONS : */
8559 /* --------------------------------- */
8560 /* There are NBRVAL values in TABLEV which stands for NBRVAL-1 intervals. */
8561 /* One searches the interval containing TPARAM by */
8562 /* dichotomy. Complexity of the algorithm : Log(n)/Log(2).(RBD). */
8564 /* ***********************************************************************
8568 /* Initialisations */
8570 /* Parameter adjustments */
8574 ibb = AdvApp2Var_SysBase::mnfndeb_();
8576 AdvApp2Var_SysBase::mgenmsg_("MMSRRE2", 7L);
8585 /* TABLEV should contain at least two values */
8592 /* TPARAM should be between extreme limits of TABLEV. */
8594 if (*tparam < tablev[1] || *tparam > tablev[*nbrval]) {
8599 /* ----------------------- SEARCH OF THE INTERVAL --------------------
8604 /* Test end of loop (found). */
8606 if (ideb + 1 == ifin) {
8611 /* Find by dichotomy on increasing values of TABLEV. */
8613 imil = (ideb + ifin) / 2;
8614 if (*tparam >= tablev[ideb] && *tparam <= tablev[imil]) {
8622 /* -------------- TEST IF TPARAM IS NOT A VALUE ---------
8623 /* ------------------------OF TABLEV UP TO EPSIL ----------------------
8627 if ((d__1 = *tparam - tablev[ideb], advapp_abs(d__1)) < *epsil) {
8631 if ((d__1 = *tparam - tablev[ifin], advapp_abs(d__1)) < *epsil) {
8636 /* --------------------------- THE END ----------------------------------
8641 AdvApp2Var_SysBase::maermsg_("MMSRRE2", iercod, 7L);
8644 AdvApp2Var_SysBase::mgsomsg_("MMSRRE2", 7L);
8649 //=======================================================================
8650 //function : mmtmave_
8652 //=======================================================================
8653 int mmtmave_(integer *nligne,
8663 /* System generated locals */
8666 /* Local variables */
8668 integer imin, imax, i__, j, k;
8673 /* ***********************************************************************
8679 /* CREATES PRODUCT G V */
8680 /* WHERE THE MATRIX IS IN FORM OF PROFILE */
8684 /* RESERVE, PRODUCT, MATRIX, PROFILE, VECTOR */
8686 /* INPUT ARGUMENTS : */
8687 /* -------------------- */
8688 /* NLIGNE : NUMBER OF LINE OF THE MATRIX */
8689 /* NCOLON : NOMBER OF COLUMN OF THE MATRIX */
8690 /* GPOSIT: TABLE OF POSITIONING OF TERMS OF STORAGE */
8691 /* GPOSIT(1,I) CONTAINS THE NUMBER of TERMS-1 ON LINE
8692 I IN THE PROFILE OF THE MATRIX */
8693 /* GPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF THE DIAGONAL TERM
8695 /* GPOSIT(3,I) CONTAINS THE INDEX COLUMN OF THE FIRST TERM OF
8696 /* PROFILE OF LINE I */
8697 /* GNSTOC : NOMBER OF TERM IN THE PROFILE OF GMATRI */
8698 /* GMATRI : MATRIX OF CONSTRAINTS IN FORM OF PROFILE */
8699 /* VECIN : INPUT VECTOR */
8701 /* OUTPUT ARGUMENTS : */
8702 /* --------------------- */
8703 /* VECOUT : VECTOR PRODUCT */
8704 /* IERCOD : ERROR CODE */
8707 /* COMMONS USED : */
8708 /* ------------------ */
8711 /* REFERENCES CALLED : */
8712 /* --------------------- */
8715 /* DESCRIPTION/NOTES/LIMITATIONS : */
8716 /* ----------------------------------- */
8718 /* ***********************************************************************
8721 /* ***********************************************************************
8726 /* ***********************************************************************
8728 /* INITIALISATIONS */
8729 /* ***********************************************************************
8732 /* Parameter adjustments */
8739 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
8741 AdvApp2Var_SysBase::mgenmsg_("MMTMAVE", 7L);
8745 /* ***********************************************************************
8748 /* ***********************************************************************
8754 for (i__ = 1; i__ <= i__1; ++i__) {
8757 for (j = 1; j <= i__2; ++j) {
8758 imin = gposit[j * 3 + 3];
8759 imax = gposit[j * 3 + 1] + gposit[j * 3 + 3] - 1;
8760 aux = gposit[j * 3 + 2] - gposit[j * 3 + 1] - imin + 1;
8761 if (imin <= i__ && i__ <= imax) {
8763 somme += gmatri[k] * vecin[j];
8766 vecout[i__] = somme;
8775 /* ***********************************************************************
8777 /* ERROR PROCESSING */
8778 /* ***********************************************************************
8782 /* ***********************************************************************
8784 /* RETURN CALLING PROGRAM */
8785 /* ***********************************************************************
8790 /* ___ DESALLOCATION, ... */
8792 AdvApp2Var_SysBase::maermsg_("MMTMAVE", iercod, 7L);
8794 AdvApp2Var_SysBase::mgsomsg_("MMTMAVE", 7L);
8799 //=======================================================================
8800 //function : mmtrpj0_
8802 //=======================================================================
8803 int mmtrpj0_(integer *ncofmx,
8813 /* System generated locals */
8814 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
8817 /* Local variables */
8819 doublereal bidon, error;
8823 /* ***********************************************************************
8828 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
8829 /* Legendre with a given precision. */
8833 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
8835 /* INPUT ARGUMENTS : */
8836 /* ------------------ */
8837 /* NCOFMX : Max Nb of coeff. of the curve (dimensioning). */
8838 /* NDIMEN : Dimension of the space. */
8839 /* NCOEFF : Degree +1 of the polynom. */
8840 /* EPSI3D : Precision required for the approximation. */
8841 /* CRVLGD : The curve the degree which of it is required to lower. */
8843 /* OUTPUT ARGUMENTS : */
8844 /* ------------------- */
8845 /* EPSTRC : Precision of the approximation. */
8846 /* NCFNEW : Degree +1 of the resulting polynom. */
8848 /* COMMONS USED : */
8849 /* ---------------- */
8851 /* REFERENCES CALLED : */
8852 /* ----------------------- */
8854 /* DESCRIPTION/NOTES/LIMITATIONS : */
8855 /* ----------------------------------- */
8857 /* ***********************************************************************
8861 /* ------- Minimum degree that can be attained : Stop at 1 (RBD) ---------
8864 /* Parameter adjustments */
8866 crvlgd_dim1 = *ncofmx;
8867 crvlgd_offset = crvlgd_dim1 + 1;
8868 crvlgd -= crvlgd_offset;
8872 /* ------------------- Init for error calculation -----------------------
8875 for (i__ = 1; i__ <= i__1; ++i__) {
8882 /* Cutting of coefficients. */
8885 /* ------ Loop on the series of Legendre :NCOEFF --> 2 (RBD) -----------
8888 for (i__ = *ncoeff; i__ >= i__1; --i__) {
8889 /* Factor of renormalization. */
8890 bidon = ((i__ - 1) * 2. + 1.) / 2.;
8891 bidon = sqrt(bidon);
8893 for (nd = 1; nd <= i__2; ++nd) {
8894 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
8898 /* Cutting is stopped if the norm becomes too great. */
8899 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
8900 if (error > *epsi3d) {
8905 /* --- Max error cumulee when the I-th coeff is removed. */
8912 /* --------------------------------- End --------------------------------
8919 //=======================================================================
8920 //function : mmtrpj2_
8922 //=======================================================================
8923 int mmtrpj2_(integer *ncofmx,
8933 /* Initialized data */
8935 static doublereal xmaxj[57] = { .9682458365518542212948163499456,
8936 .986013297183269340427888048593603,
8937 1.07810420343739860362585159028115,
8938 1.17325804490920057010925920756025,
8939 1.26476561266905634732910520370741,
8940 1.35169950227289626684434056681946,
8941 1.43424378958284137759129885012494,
8942 1.51281316274895465689402798226634,
8943 1.5878364329591908800533936587012,
8944 1.65970112228228167018443636171226,
8945 1.72874345388622461848433443013543,
8946 1.7952515611463877544077632304216,
8947 1.85947199025328260370244491818047,
8948 1.92161634324190018916351663207101,
8949 1.98186713586472025397859895825157,
8950 2.04038269834980146276967984252188,
8951 2.09730119173852573441223706382076,
8952 2.15274387655763462685970799663412,
8953 2.20681777186342079455059961912859,
8954 2.25961782459354604684402726624239,
8955 2.31122868752403808176824020121524,
8956 2.36172618435386566570998793688131,
8957 2.41117852396114589446497298177554,
8958 2.45964731268663657873849811095449,
8959 2.50718840313973523778244737914028,
8960 2.55385260994795361951813645784034,
8961 2.59968631659221867834697883938297,
8962 2.64473199258285846332860663371298,
8963 2.68902863641518586789566216064557,
8964 2.73261215675199397407027673053895,
8965 2.77551570192374483822124304745691,
8966 2.8177699459714315371037628127545,
8967 2.85940333797200948896046563785957,
8968 2.90044232019793636101516293333324,
8969 2.94091151970640874812265419871976,
8970 2.98083391718088702956696303389061,
8971 3.02023099621926980436221568258656,
8972 3.05912287574998661724731962377847,
8973 3.09752842783622025614245706196447,
8974 3.13546538278134559341444834866301,
8975 3.17295042316122606504398054547289,
8976 3.2099992681699613513775259670214,
8977 3.24662674946606137764916854570219,
8978 3.28284687953866689817670991319787,
8979 3.31867291347259485044591136879087,
8980 3.35411740487202127264475726990106,
8981 3.38919225660177218727305224515862,
8982 3.42390876691942143189170489271753,
8983 3.45827767149820230182596660024454,
8984 3.49230918177808483937957161007792,
8985 3.5260130200285724149540352829756,
8986 3.55939845146044235497103883695448,
8987 3.59247431368364585025958062194665,
8988 3.62524904377393592090180712976368,
8989 3.65773070318071087226169680450936,
8990 3.68992700068237648299565823810245,
8991 3.72184531357268220291630708234186 };
8993 /* System generated locals */
8994 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
8997 /* Local variables */
8999 doublereal bidon, error;
9001 doublereal bid, eps1;
9004 /* ***********************************************************************
9009 /* Lower the degree of a curve defined on (-1,1) in the direction of */
9010 /* Legendre with a given precision. */
9014 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
9016 /* INPUT ARGUMENTS : */
9017 /* ------------------ */
9018 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9019 /* NDIMEN : Dimension of the space. */
9020 /* NCOEFF : Degree +1 of the polynom. */
9021 /* EPSI3D : Precision required for the approximation. */
9022 /* CRVLGD : The curve the degree which of will be lowered. */
9024 /* OUTPUT ARGUMENTS : */
9025 /* ------------------- */
9026 /* YCVMAX : Auxiliary table (error max on each dimension).
9028 /* EPSTRC : Precision of the approximation. */
9029 /* NCFNEW : Degree +1 of the resulting polynom. */
9031 /* COMMONS USED : */
9032 /* ---------------- */
9034 /* REFERENCES CALLED : */
9035 /* ----------------------- */
9037 /* DESCRIPTION/NOTES/LIMITATIONS : */
9038 /* ----------------------------------- */
9040 /* ***********************************************************************
9044 /* Parameter adjustments */
9046 crvlgd_dim1 = *ncofmx;
9047 crvlgd_offset = crvlgd_dim1 + 1;
9048 crvlgd -= crvlgd_offset;
9054 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9058 /* Init for calculation of error. */
9060 for (i__ = 1; i__ <= i__1; ++i__) {
9067 /* Cutting of coefficients. */
9070 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9073 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9074 /* Factor of renormalization. */
9075 bidon = xmaxj[i__ - ncut];
9077 for (nd = 1; nd <= i__2; ++nd) {
9078 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
9082 /* One stops to cut if the norm becomes too great. */
9083 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9084 if (error > *epsi3d) {
9089 /* --- Max error cumulated when the I-th coeff is removed. */
9096 /* ------- Cutting of zero coeffs of interpolation (RBD) -------
9100 if (*ncfnew == ia) {
9101 AdvApp2Var_MathBase::mmeps1_(&eps1);
9102 for (i__ = ia; i__ >= 2; --i__) {
9105 for (nd = 1; nd <= i__1; ++nd) {
9106 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
9115 /* --- If all coeffs can be removed, this is a point. */
9119 /* --------------------------------- End --------------------------------
9126 //=======================================================================
9127 //function : mmtrpj4_
9129 //=======================================================================
9130 int mmtrpj4_(integer *ncofmx,
9139 /* Initialized data */
9141 static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
9142 1.05299572648705464724876659688996,
9143 1.0949715351434178709281698645813,
9144 1.15078388379719068145021100764647,
9145 1.2094863084718701596278219811869,
9146 1.26806623151369531323304177532868,
9147 1.32549784426476978866302826176202,
9148 1.38142537365039019558329304432581,
9149 1.43575531950773585146867625840552,
9150 1.48850442653629641402403231015299,
9151 1.53973611681876234549146350844736,
9152 1.58953193485272191557448229046492,
9153 1.63797820416306624705258190017418,
9154 1.68515974143594899185621942934906,
9155 1.73115699602477936547107755854868,
9156 1.77604489805513552087086912113251,
9157 1.81989256661534438347398400420601,
9158 1.86276344480103110090865609776681,
9159 1.90471563564740808542244678597105,
9160 1.94580231994751044968731427898046,
9161 1.98607219357764450634552790950067,
9162 2.02556989246317857340333585562678,
9163 2.06433638992049685189059517340452,
9164 2.10240936014742726236706004607473,
9165 2.13982350649113222745523925190532,
9166 2.17661085564771614285379929798896,
9167 2.21280102016879766322589373557048,
9168 2.2484214321456956597803794333791,
9169 2.28349755104077956674135810027654,
9170 2.31805304852593774867640120860446,
9171 2.35210997297725685169643559615022,
9172 2.38568889602346315560143377261814,
9173 2.41880904328694215730192284109322,
9174 2.45148841120796359750021227795539,
9175 2.48374387161372199992570528025315,
9176 2.5155912654873773953959098501893,
9177 2.54704548720896557684101746505398,
9178 2.57812056037881628390134077704127,
9179 2.60882970619319538196517982945269,
9180 2.63918540521920497868347679257107,
9181 2.66919945330942891495458446613851,
9182 2.69888301230439621709803756505788,
9183 2.72824665609081486737132853370048,
9184 2.75730041251405791603760003778285,
9185 2.78605380158311346185098508516203,
9186 2.81451587035387403267676338931454,
9187 2.84269522483114290814009184272637,
9188 2.87060005919012917988363332454033,
9189 2.89823818258367657739520912946934,
9190 2.92561704377132528239806135133273,
9191 2.95274375377994262301217318010209,
9192 2.97962510678256471794289060402033,
9193 3.00626759936182712291041810228171,
9194 3.03267744830655121818899164295959,
9195 3.05886060707437081434964933864149 };
9197 /* System generated locals */
9198 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
9201 /* Local variables */
9203 doublereal bidon, error;
9205 doublereal bid, eps1;
9209 /* ***********************************************************************
9214 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
9215 /* Legendre with a given precision. */
9219 /* LEGENDRE, POLYGON, TRONCATION, CURVE, SMOOTHING. */
9221 /* INPUT ARGUMENTS : */
9222 /* ------------------ */
9223 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9224 /* NDIMEN : Dimension of the space. */
9225 /* NCOEFF : Degree +1 of the polynom. */
9226 /* EPSI3D : Precision required for the approximation. */
9227 /* CRVLGD : The curve which wishes to lower the degree. */
9229 /* OUTPUT ARGUMENTS : */
9230 /* ------------------- */
9231 /* YCVMAX : Auxiliary table (max error on each dimension).
9233 /* EPSTRC : Precision of the approximation. */
9234 /* NCFNEW : Degree +1 of the resulting polynom. */
9236 /* COMMONS USED : */
9237 /* ---------------- */
9239 /* REFERENCES CALLED : */
9240 /* ----------------------- */
9242 /* DESCRIPTION/NOTES/LIMITATIONS : */
9243 /* ----------------------------------- */
9245 /* ***********************************************************************
9249 /* Parameter adjustments */
9251 crvlgd_dim1 = *ncofmx;
9252 crvlgd_offset = crvlgd_dim1 + 1;
9253 crvlgd -= crvlgd_offset;
9259 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9263 /* Init for error calculation. */
9265 for (i__ = 1; i__ <= i__1; ++i__) {
9272 /* Cutting of coefficients. */
9275 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9278 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9279 /* Factor of renormalization. */
9280 bidon = xmaxj[i__ - ncut];
9282 for (nd = 1; nd <= i__2; ++nd) {
9283 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
9287 /* Stop cutting if the norm becomes too great. */
9288 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9289 if (error > *epsi3d) {
9294 /* -- Error max cumulated when the I-eme coeff is removed. */
9301 /* ------- Cutting of zero coeffs of the pole of interpolation (RBD) -------
9305 if (*ncfnew == ia) {
9306 AdvApp2Var_MathBase::mmeps1_(&eps1);
9307 for (i__ = ia; i__ >= 2; --i__) {
9310 for (nd = 1; nd <= i__1; ++nd) {
9311 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
9320 /* --- If all coeffs can be removed, this is a point. */
9324 /* --------------------------------- End --------------------------------
9331 //=======================================================================
9332 //function : mmtrpj6_
9334 //=======================================================================
9335 int mmtrpj6_(integer *ncofmx,
9345 /* Initialized data */
9347 static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
9348 1.11626917091567929907256116528817,
9349 1.1327140810290884106278510474203,
9350 1.1679452722668028753522098022171,
9351 1.20910611986279066645602153641334,
9352 1.25228283758701572089625983127043,
9353 1.29591971597287895911380446311508,
9354 1.3393138157481884258308028584917,
9355 1.3821288728999671920677617491385,
9356 1.42420414683357356104823573391816,
9357 1.46546895108549501306970087318319,
9358 1.50590085198398789708599726315869,
9359 1.54550385142820987194251585145013,
9360 1.58429644271680300005206185490937,
9361 1.62230484071440103826322971668038,
9362 1.65955905239130512405565733793667,
9363 1.69609056468292429853775667485212,
9364 1.73193098017228915881592458573809,
9365 1.7671112206990325429863426635397,
9366 1.80166107681586964987277458875667,
9367 1.83560897003644959204940535551721,
9368 1.86898184653271388435058371983316,
9369 1.90180515174518670797686768515502,
9370 1.93410285411785808749237200054739,
9371 1.96589749778987993293150856865539,
9372 1.99721027139062501070081653790635,
9373 2.02806108474738744005306947877164,
9374 2.05846864831762572089033752595401,
9375 2.08845055210580131460156962214748,
9376 2.11802334209486194329576724042253,
9377 2.14720259305166593214642386780469,
9378 2.17600297710595096918495785742803,
9379 2.20443832785205516555772788192013,
9380 2.2325216999457379530416998244706,
9381 2.2602654243075083168599953074345,
9382 2.28768115912702794202525264301585,
9383 2.3147799369092684021274946755348,
9384 2.34157220782483457076721300512406,
9385 2.36806787963276257263034969490066,
9386 2.39427635443992520016789041085844,
9387 2.42020656255081863955040620243062,
9388 2.44586699364757383088888037359254,
9389 2.47126572552427660024678584642791,
9390 2.49641045058324178349347438430311,
9391 2.52130850028451113942299097584818,
9392 2.54596686772399937214920135190177,
9393 2.5703922285006754089328998222275,
9394 2.59459096001908861492582631591134,
9395 2.61856915936049852435394597597773,
9396 2.64233265984385295286445444361827,
9397 2.66588704638685848486056711408168,
9398 2.68923766976735295746679957665724,
9399 2.71238965987606292679677228666411 };
9401 /* System generated locals */
9402 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
9405 /* Local variables */
9407 doublereal bidon, error;
9409 doublereal bid, eps1;
9413 /* ***********************************************************************
9418 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
9419 /* Legendre to a given precision. */
9423 /* LEGENDRE,POLYGON,TRUNCATION,CURVE,SMOOTHING. */
9425 /* INPUT ARGUMENTS : */
9426 /* ------------------ */
9427 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9428 /* NDIMEN : Dimension of the space. */
9429 /* NCOEFF : Degree +1 of the polynom. */
9430 /* EPSI3D : Precision required for the approximation. */
9431 /* CRVLGD : The curve the degree which of will be lowered. */
9433 /* OUTPUT ARGUMENTS : */
9434 /* ------------------- */
9435 /* YCVMAX : Auxiliary table (max error on each dimension).
9436 /* EPSTRC : Precision of the approximation. */
9437 /* NCFNEW : Degree +1 of the resulting polynom. */
9439 /* COMMONS USED : */
9440 /* ---------------- */
9442 /* REFERENCES CALLED : */
9443 /* ----------------------- */
9445 /* DESCRIPTION/NOTES/LIMITATIONS : */
9446 /* ----------------------------------- */
9448 /* ***********************************************************************
9452 /* Parameter adjustments */
9454 crvlgd_dim1 = *ncofmx;
9455 crvlgd_offset = crvlgd_dim1 + 1;
9456 crvlgd -= crvlgd_offset;
9462 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9466 /* Init for error calculation. */
9468 for (i__ = 1; i__ <= i__1; ++i__) {
9475 /* Cutting of coefficients. */
9478 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9481 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9482 /* Factor of renormalization. */
9483 bidon = xmaxj[i__ - ncut];
9485 for (nd = 1; nd <= i__2; ++nd) {
9486 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
9490 /* Stop cutting if the norm becomes too great. */
9491 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9492 if (error > *epsi3d) {
9497 /* --- Max error cumulated when the I-th coeff is removed. */
9504 /* ------- Cutting of zero coeff. of the pole of interpolation (RBD) -------
9508 if (*ncfnew == ia) {
9509 AdvApp2Var_MathBase::mmeps1_(&eps1);
9510 for (i__ = ia; i__ >= 2; --i__) {
9513 for (nd = 1; nd <= i__1; ++nd) {
9514 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
9523 /* --- If all coeffs can be removed, this is a point. */
9527 /* --------------------------------- End --------------------------------
9534 //=======================================================================
9535 //function : AdvApp2Var_MathBase::mmtrpjj_
9537 //=======================================================================
9538 int AdvApp2Var_MathBase::mmtrpjj_(integer *ncofmx,
9548 /* System generated locals */
9549 integer crvlgd_dim1, crvlgd_offset;
9551 /* Local variables */
9555 /* ***********************************************************************
9560 /* Lower the degree of a curve defined on (-1,1) in the direction of */
9561 /* Legendre with a given precision. */
9565 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
9567 /* INPUT ARGUMENTS : */
9568 /* ------------------ */
9569 /* NCOFMX : Max Nb coeff. of the curve (dimensioning). */
9570 /* NDIMEN : Dimension of the space. */
9571 /* NCOEFF : Degree +1 of the polynom. */
9572 /* EPSI3D : Precision required for the approximation. */
9573 /* IORDRE : Order of continuity at the extremities. */
9574 /* CRVLGD : The curve the degree which of should be lowered. */
9576 /* OUTPUT ARGUMENTS : */
9577 /* ------------------- */
9578 /* ERRMAX : Precision of the approximation. */
9579 /* NCFNEW : Degree +1 of the resulting polynom. */
9581 /* COMMONS USED : */
9582 /* ---------------- */
9584 /* REFERENCES CALLED : */
9585 /* ----------------------- */
9587 /* DESCRIPTION/NOTES/LIMITATIONS : */
9588 /* ----------------------------------- */
9590 /* ***********************************************************************
9594 /* Parameter adjustments */
9596 crvlgd_dim1 = *ncofmx;
9597 crvlgd_offset = crvlgd_dim1 + 1;
9598 crvlgd -= crvlgd_offset;
9601 ia = (*iordre + 1) << 1;
9604 mmtrpj0_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9605 ycvmax[1], errmax, ncfnew);
9606 } else if (ia == 2) {
9607 mmtrpj2_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9608 ycvmax[1], errmax, ncfnew);
9609 } else if (ia == 4) {
9610 mmtrpj4_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9611 ycvmax[1], errmax, ncfnew);
9613 mmtrpj6_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9614 ycvmax[1], errmax, ncfnew);
9617 /* ------------------------ End -----------------------------------------
9623 //=======================================================================
9624 //function : AdvApp2Var_MathBase::mmunivt_
9626 //=======================================================================
9627 int AdvApp2Var_MathBase::mmunivt_(integer *ndimen,
9634 doublereal c_b2 = 10.;
9636 /* System generated locals */
9640 /* Local variables */
9641 integer nchif, iunit, izero;
9650 /* ***********************************************************************
9655 /* CALCULATE THE NORMAL VECTOR BASING ON ANY VECTOR */
9656 /* WITH PRECISION GIVEN BY THE USER. */
9660 /* ALL, MATH_ACCES :: */
9661 /* VECTEUR&, NORMALISATION, &VECTEUR */
9663 /* INPUT ARGUMENTS : */
9664 /* ------------------ */
9665 /* NDIMEN : DIMENSION OF THE SPACE */
9666 /* VECTOR : VECTOR TO BE NORMED */
9667 /* EPSILN : EPSILON BELOW WHICH IT IS CONSIDERED THAT THE */
9668 /* NORM OF THE VECTOR IS NULL. IF EPSILN<=0, A DEFAULT VALUE */
9669 /* IS IMPOSED (10.D-17 ON VAX). */
9671 /* OUTPUT ARGUMENTS : */
9672 /* ------------------- */
9673 /* VECNRM : NORMED VECTOR */
9674 /* IERCOD 101 : THE VECTOR IS NULL UP TO EPSILN. */
9677 /* COMMONS USED : */
9678 /* ---------------- */
9680 /* REFERENCES CALLED : */
9681 /* ----------------------- */
9683 /* DESCRIPTION/NOTES/LIMITATIONS : */
9684 /* ----------------------------------- */
9685 /* VECTOR and VECNRM can be identic. */
9687 /* The norm of vector is calculated and each component is divided by
9688 /* this norm. After this it is checked if all componentes of the */
9689 /* vector except for one cost 0 with machine precision. In */
9690 /* this case the quasi-null components are set to 0.D0. */
9692 /* ***********************************************************************
9696 /* Parameter adjustments */
9703 /* -------- Precision by default : zero machine 10.D-17 on Vax ------
9706 AdvApp2Var_SysBase::maovsr8_(&nchif);
9707 if (*epsiln <= 0.) {
9709 eps0 = AdvApp2Var_MathBase::pow__di(&c_b2, &i__1);
9714 /* ------------------------- Calculation of the norm --------------------
9717 vnorm = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vector[1]);
9718 if (vnorm <= eps0) {
9719 AdvApp2Var_SysBase::mvriraz_(ndimen, &vecnrm[1]);
9724 /* ---------------------- Calculation of the vector norm ---------------
9728 i__1 = (-nchif - 1) / 2;
9729 eps0 = AdvApp2Var_MathBase::pow__di(&c_b2, &i__1);
9731 for (ii = 1; ii <= i__1; ++ii) {
9732 vecnrm[ii] = vector[ii] / vnorm;
9733 if ((d__1 = vecnrm[ii], advapp_abs(d__1)) <= eps0) {
9741 /* ------ Case when all coordinates except for one are almost null ----
9743 /* ------------- then one of coordinates costs 1.D0 or -1.D0 --------
9746 if (izero == *ndimen - 1) {
9747 bid = vecnrm[iunit];
9749 for (ii = 1; ii <= i__1; ++ii) {
9756 vecnrm[iunit] = -1.;
9760 /* -------------------------------- The end -----------------------------
9767 //=======================================================================
9768 //function : AdvApp2Var_MathBase::mmveps3_
9770 //=======================================================================
9771 int AdvApp2Var_MathBase::mmveps3_(doublereal *eps03)
9773 /* Initialized data */
9775 static char nomprg[8+1] = "MMEPS1 ";
9781 /************************************************************************
9786 /* Extraction of EPS1 from COMMON MPRCSN. */
9790 /* MPRCSN,PRECISON,EPS3. */
9792 /* INPUT ARGUMENTS : */
9793 /* ------------------ */
9796 /* OUTPUT ARGUMENTS : */
9797 /* ------------------- */
9798 /* EPS3 : space zero of the denominator (10**-9) */
9799 /* EPS3 should value 10**-15 */
9801 /* COMMONS USED : */
9802 /* ---------------- */
9804 /* REFERENCES CALLED : */
9805 /* ----------------------- */
9807 /* DESCRIPTION/NOTES/LIMITATIONS : */
9808 /* ----------------------------------- */
9811 /* ***********************************************************************
9816 /* ***********************************************************************
9821 /* GIVES TOLERANCES OF NULLITY IN STRIM */
9822 /* AND LIMITS OF ITERATIVE PROCESSES */
9824 /* GENERAL CONTEXT, MODIFIABLE BY THE UTILISER */
9828 /* PARAMETER , TOLERANCE */
9830 /* DESCRIPTION/NOTES/LIMITATIONS : */
9831 /* ----------------------------------- */
9832 /* INITIALISATION : PROFILE , **VIA MPRFTX** AT INPUT IN STRIM*/
9833 /* LOADING OF DEFAULT VALUES OF THE PROFILE IN MPRFTX AT INPUT*/
9834 /* IN STRIM. THEY ARE PRESERVED IN THE LOCAL VARIABLES OF MPRFTX */
9836 /* RESET DEFAULT VALUES : MDFINT */
9837 /* MODIFICATION INTERACTIVE BY THE USER : MDBINT */
9839 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
9840 /* MEPSPB ... EPS3,EPS4 */
9841 /* MEPSLN ... EPS2, NITERM , NITERR */
9842 /* MEPSNR ... EPS2 , NITERM */
9843 /* MITERR ... NITERR */
9846 /* ***********************************************************************
9849 /* NITERM : MAX NB OF ITERATIONS */
9850 /* NITERR : NB OF RAPID ITERATIONS */
9851 /* EPS1 : TOLERANCE OF 3D NULL DISTANCE */
9852 /* EPS2 : TOLERANCE OF ZERO PARAMETRIC DISTANCE */
9853 /* EPS3 : TOLERANCE TO AVOID DIVISION BY 0.. */
9854 /* EPS4 : TOLERANCE ANGULAR */
9858 /* ***********************************************************************
9861 ibb = AdvApp2Var_SysBase::mnfndeb_();
9863 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
9866 *eps03 = mmprcsn_.eps3;
9871 //=======================================================================
9872 //function : AdvApp2Var_MathBase::mmvncol_
9874 //=======================================================================
9875 int AdvApp2Var_MathBase::mmvncol_(integer *ndimen,
9881 /* System generated locals */
9884 /* Local variables */
9887 doublereal vaux1[3], vaux2[3];
9893 /* ***********************************************************************
9898 /* CALCULATE A VECTOR NON-COLINEAR TO A GIVEN NON-NULL VECTOR */
9902 /* PUBLIC, VECTOR, FREE */
9904 /* INPUT ARGUMENTS : */
9905 /* -------------------- */
9906 /* ndimen : dimension of the space */
9907 /* vecin : input vector */
9909 /* OUTPUT ARGUMENTS : */
9910 /* --------------------- */
9912 /* vecout : vector non colinear to vecin */
9914 /* COMMONS USED : */
9915 /* ------------------ */
9918 /* REFERENCES CALLED : */
9919 /* --------------------- */
9922 /* DESCRIPTION/NOTES/LIMITATIONS : */
9923 /* ----------------------------------- */
9925 /* ***********************************************************************
9928 /* ***********************************************************************
9933 /* ***********************************************************************
9935 /* INITIALISATIONS */
9936 /* ***********************************************************************
9939 /* Parameter adjustments */
9944 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
9946 AdvApp2Var_SysBase::mgenmsg_("MMVNCOL", 7L);
9950 /* ***********************************************************************
9953 /* ***********************************************************************
9956 if (*ndimen <= 1 || *ndimen > 3) {
9962 while(d__ <= *ndimen) {
9963 if (vecin[d__] == 0.) {
9968 if (aux == *ndimen) {
9973 for (d__ = 1; d__ <= 3; ++d__) {
9974 vaux1[d__ - 1] = 0.;
9977 for (d__ = 1; d__ <= i__1; ++d__) {
9978 vaux1[d__ - 1] = vecin[d__];
9979 vaux2[d__ - 1] = vecin[d__];
9988 vaux2[d__ - 1] += 1;
9989 valaux = vaux1[1] * vaux2[2] - vaux1[2] * vaux2[1];
9991 valaux = vaux1[2] * vaux2[0] - vaux1[0] * vaux2[2];
9993 valaux = vaux1[0] * vaux2[1] - vaux1[1] * vaux2[0];
10008 for (d__ = 1; d__ <= i__1; ++d__) {
10009 vecout[d__] = vaux2[d__ - 1];
10014 /* ***********************************************************************
10016 /* ERROR PROCESSING */
10017 /* ***********************************************************************
10026 /* ***********************************************************************
10028 /* RETURN CALLING PROGRAM */
10029 /* ***********************************************************************
10035 AdvApp2Var_SysBase::maermsg_("MMVNCOL", iercod, 7L);
10037 AdvApp2Var_SysBase::mgsomsg_("MMVNCOL", 7L);
10042 //=======================================================================
10043 //function : AdvApp2Var_MathBase::mmwprcs_
10045 //=======================================================================
10046 void AdvApp2Var_MathBase::mmwprcs_(doublereal *epsil1,
10047 doublereal *epsil2,
10048 doublereal *epsil3,
10049 doublereal *epsil4,
10056 /* ***********************************************************************
10061 /* ACCESS IN WRITING FOR COMMON MPRCSN */
10067 /* INPUT ARGUMENTS : */
10068 /* -------------------- */
10069 /* EPSIL1 : TOLERANCE OF 3D NULL DISTANCE */
10070 /* EPSIL2 : TOLERANCE OF PARAMETRIC NULL DISTANCE */
10071 /* EPSIL3 : TOLERANCE TO AVOID DIVISION BY 0.. */
10072 /* EPSIL4 : ANGULAR TOLERANCE */
10073 /* NITER1 : MAX NB OF ITERATIONS */
10074 /* NITER2 : NB OF RAPID ITERATIONS */
10076 /* OUTPUT ARGUMENTS : */
10077 /* --------------------- */
10080 /* COMMONS USED : */
10081 /* ------------------ */
10084 /* REFERENCES CALLED : */
10085 /* --------------------- */
10088 /* DESCRIPTION/NOTES/LIMITATIONS : */
10089 /* ----------------------------------- */
10092 /* ***********************************************************************
10095 /* ***********************************************************************
10099 /* ***********************************************************************
10101 /* INITIALIZATIONS */
10102 /* ***********************************************************************
10105 /* ***********************************************************************
10108 /* ***********************************************************************
10111 /* ***********************************************************************
10116 /* GIVES TOLERANCES OF NULLITY IN STRIM */
10117 /* AND LIMITS OF ITERATIVE PROCESSES */
10119 /* GENERAL CONTEXT, MODIFIABLE BY THE UTILISER */
10123 /* PARAMETER , TOLERANCE */
10125 /* DESCRIPTION/NOTES/LIMITATIONS : */
10126 /* ----------------------------------- */
10127 /* INITIALISATION : PROFILE , **VIA MPRFTX** AT INPUT IN STRIM*/
10128 /* LOADING OF DEFAULT VALUES OF THE PROFILE IN MPRFTX AT INPUT*/
10129 /* IN STRIM. THEY ARE PRESERVED IN THE LOCAL VARIABLES OF MPRFTX */
10131 /* RESET DEFAULT VALUES : MDFINT */
10132 /* MODIFICATION INTERACTIVE BY THE USER : MDBINT */
10134 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
10135 /* MEPSPB ... EPS3,EPS4 */
10136 /* MEPSLN ... EPS2, NITERM , NITERR */
10137 /* MEPSNR ... EPS2 , NITERM */
10138 /* MITERR ... NITERR */
10141 /* ***********************************************************************
10144 /* NITERM : MAX NB OF ITERATIONS */
10145 /* NITERR : NB OF RAPID ITERATIONS */
10146 /* EPS1 : TOLERANCE OF 3D NULL DISTANCE */
10147 /* EPS2 : TOLERANCE OF ZERO PARAMETRIC DISTANCE */
10148 /* EPS3 : TOLERANCE TO AVOID DIVISION BY 0.. */
10149 /* EPS4 : TOLERANCE ANGULAR */
10152 /* ***********************************************************************
10154 mmprcsn_.eps1 = *epsil1;
10155 mmprcsn_.eps2 = *epsil2;
10156 mmprcsn_.eps3 = *epsil3;
10157 mmprcsn_.eps4 = *epsil4;
10158 mmprcsn_.niterm = *niter1;
10159 mmprcsn_.niterr = *niter2;
10164 //=======================================================================
10165 //function : AdvApp2Var_MathBase::pow__di
10167 //=======================================================================
10168 doublereal AdvApp2Var_MathBase::pow__di (doublereal *x,
10172 register integer ii ;
10173 doublereal result ;
10176 if ( *n > 0 ) {absolute = *n;}
10177 else {absolute = -*n;}
10178 /* System generated locals */
10179 for(ii = 0 ; ii < absolute ; ii++) {
10183 result = 1.0e0 / result ;
10189 /* **********************************************************************
10194 /* Calculate integer function power not obligatory in the most efficient way ;
10201 /* INPUT ARGUMENTS : */
10202 /* ------------------ */
10203 /* X : argument of X**N */
10206 /* OUTPUT ARGUMENTS : */
10207 /* ------------------- */
10210 /* COMMONS USED : */
10211 /* ---------------- */
10213 /* REFERENCES CALLED : */
10214 /* ----------------------- */
10216 /* DESCRIPTION/NOTES/LIMITATIONS : */
10217 /* ----------------------------------- */
10220 /* ***********************************************************************/
10222 //=======================================================================
10223 //function : pow__ii
10225 //=======================================================================
10226 integer pow__ii(integer *x,
10230 register integer ii ;
10234 if ( *n > 0 ) {absolute = *n;}
10235 else {absolute = -*n;}
10236 /* System generated locals */
10237 for(ii = 0 ; ii < absolute ; ii++) {
10241 result = 1 / result ;
10247 /* **********************************************************************
10249 /* **********************************************************************
10254 /* Calculate integer function power not obligatory in the most efficient way ;
10261 /* INPUT ARGUMENTS : */
10262 /* ------------------ */
10263 /* X : argument of X**N */
10266 /* OUTPUT ARGUMENTS : */
10267 /* ------------------- */
10270 /* COMMONS USED : */
10271 /* ---------------- */
10273 /* REFERENCES CALLED : */
10274 /* ----------------------- */
10276 /* DESCRIPTION/NOTES/LIMITATIONS : */
10277 /* ----------------------------------- */
10280 /* ***********************************************************************/
10282 //=======================================================================
10283 //function : AdvApp2Var_MathBase::msc_
10285 //=======================================================================
10286 doublereal AdvApp2Var_MathBase::msc_(integer *ndimen,
10287 doublereal *vecte1,
10288 doublereal *vecte2)
10291 /* System generated locals */
10293 doublereal ret_val;
10295 /* Local variables */
10301 /************************************************************************
10306 /* Calculate the scalar product of 2 vectors in the space */
10307 /* of dimension NDIMEN. */
10311 /* PRODUCT MSCALAIRE. */
10313 /* INPUT ARGUMENTS : */
10314 /* ------------------ */
10315 /* NDIMEN : Dimension of the space. */
10316 /* VECTE1,VECTE2: Vectors. */
10318 /* OUTPUT ARGUMENTS : */
10319 /* ------------------- */
10321 /* COMMONS USED : */
10322 /* ---------------- */
10324 /* REFERENCES CALLED : */
10325 /* ----------------------- */
10327 /* DESCRIPTION/NOTES/LIMITATIONS : */
10328 /* ----------------------------------- */
10331 /* ***********************************************************************
10335 /* PRODUIT MSCALAIRE */
10336 /* Parameter adjustments */
10340 /* Function Body */
10344 for (i__ = 1; i__ <= i__1; ++i__) {
10345 x += vecte1[i__] * vecte2[i__];
10350 /* ----------------------------------- THE END --------------------------
10356 //=======================================================================
10357 //function : mvcvin2_
10359 //=======================================================================
10360 int mvcvin2_(integer *ncoeff,
10361 doublereal *crvold,
10362 doublereal *crvnew,
10366 /* System generated locals */
10367 integer i__1, i__2;
10369 /* Local variables */
10370 integer m1jm1, ncfm1, j, k;
10372 doublereal cij1, cij2;
10376 /************************************************************************
10381 /* INVERSION OF THE PARAMETERS ON CURVE 2D. */
10385 /* CURVE,2D,INVERSION,PARAMETER. */
10387 /* INPUT ARGUMENTS : */
10388 /* ------------------ */
10389 /* NCOEFF : NB OF COEFF OF THE CURVE. */
10390 /* CRVOLD : CURVE OF ORIGIN */
10392 /* OUTPUT ARGUMENTS : */
10393 /* ------------------- */
10394 /* CRVNEW : THE RESULTING CURVE AFTER CHANGE OF T BY 1-T */
10395 /* IERCOD : 0 OK, */
10396 /* 10 NB OF COEFF NULL OR TOO GREAT. */
10398 /* COMMONS USED : */
10399 /* ---------------- */
10402 /* REFERENCES CALLED : */
10403 /* ---------------------- */
10405 /* DESCRIPTION/NOTES/LIMITATIONS : */
10406 /* ----------------------------------- */
10407 /* THE FOLLOWING CALL IS ABSOLUTELY LEGAL : */
10408 /* CALL MVCVIN2(NCOEFF,CURVE,CURVE,IERCOD), THE TABLE CURVE */
10409 /* BECOMES INPUT AND OUTPUT ARGUMENT (RBD). */
10410 /* BECAUSE OF MCCNP, THE NB OF COEFF OF THE CURVE IS LIMITED TO */
10411 /* NDGCNP+1 = 61. */
10414 /* ***********************************************************************
10418 /* **********************************************************************
10423 /* Serves to provide coefficients of the binome (triangle of Pascal). */
10427 /* Coeff of binome from 0 to 60. read only . init par block data */
10429 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
10430 /* ----------------------------------- */
10431 /* The coefficients of the binome form a triangular matrix. */
10432 /* This matrix is completed in table CNP by transposition. */
10433 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
10435 /* Initialization is done by block-data MMLLL09.RES, */
10436 /* created by program MQINICNP.FOR (see the team (AC) ). */
10440 /* **********************************************************************
10445 /* ***********************************************************************
10448 /* Parameter adjustments */
10452 /* Function Body */
10453 if (*ncoeff < 1 || *ncoeff - 1 > 60) {
10460 /* CONSTANT TERM OF THE NEW CURVE */
10465 for (k = 2; k <= i__1; ++k) {
10466 cij1 += crvold[(k << 1) + 1];
10467 cij2 += crvold[(k << 1) + 2];
10471 if (*ncoeff == 1) {
10475 /* INTERMEDIARY POWERS OF THE PARAMETER */
10477 ncfm1 = *ncoeff - 1;
10480 for (j = 2; j <= i__1; ++j) {
10482 cij1 = crvold[(j << 1) + 1];
10483 cij2 = crvold[(j << 1) + 2];
10485 for (k = j + 1; k <= i__2; ++k) {
10486 bid = mmcmcnp_.cnp[k - 1 + (j - 1) * 61];
10487 cij1 += crvold[(k << 1) + 1] * bid;
10488 cij2 += crvold[(k << 1) + 2] * bid;
10490 crvnew[(j << 1) + 1] = cij1 * m1jm1;
10491 crvnew[(j << 1) + 2] = cij2 * m1jm1;
10494 /* TERM OF THE HIGHEST DEGREE */
10496 crvnew[(*ncoeff << 1) + 1] = -crvold[(*ncoeff << 1) + 1] * m1jm1;
10497 crvnew[(*ncoeff << 1) + 2] = -crvold[(*ncoeff << 1) + 2] * m1jm1;
10501 AdvApp2Var_SysBase::maermsg_("MVCVIN2", iercod, 7L);
10506 //=======================================================================
10507 //function : mvcvinv_
10509 //=======================================================================
10510 int mvcvinv_(integer *ncoeff,
10511 doublereal *crvold,
10512 doublereal *crvnew,
10516 /* System generated locals */
10517 integer i__1, i__2;
10519 /* Local variables */
10520 integer m1jm1, ncfm1, j, k;
10522 //extern /* Subroutine */ int maermsg_();
10523 doublereal cij1, cij2, cij3;
10526 /* **********************************************************************
10531 /* INVERSION OF THE PARAMETER ON A CURBE 3D (I.E. INVERSION */
10532 /* OF THE DIRECTION OF PARSING). */
10536 /* CURVE,INVERSION,PARAMETER. */
10538 /* INPUT ARGUMENTS : */
10539 /* ------------------ */
10540 /* NCOEFF : NB OF COEFF OF THE CURVE. */
10541 /* CRVOLD : CURVE OF ORIGIN */
10543 /* OUTPUT ARGUMENTS : */
10544 /* ------------------- */
10545 /* CRVNEW : RESULTING CURVE AFTER CHANGE OF T INTO 1-T */
10546 /* IERCOD : 0 OK, */
10547 /* 10 NB OF COEFF NULL OR TOO GREAT. */
10549 /* COMMONS USED : */
10550 /* ---------------- */
10553 /* REFERENCES CALLED : */
10554 /* ---------------------- */
10556 /* DESCRIPTION/NOTES/LIMITATIONS : */
10557 /* ----------------------------------- */
10558 /* THE FOLLOWING CALL IS ABSOLUTELY LEGAL : */
10559 /* CALL MVCVINV(NCOEFF,CURVE,CURVE,IERCOD), TABLE CURVE */
10560 /* BECOMES INPUT AND OUTPUT ARGUMENT (RBD). */
10561 /* THE NUMBER OF COEFF OF THE CURVE IS LIMITED TO NDGCNP+1 = 61 */
10562 /* BECAUSE OF USE OF COMMON MCCNP. */
10564 /* ***********************************************************************
10567 /* **********************************************************************
10572 /* Serves to provide the binomial coefficients (triangle of Pascal). */
10576 /* Binomial Coeff from 0 to 60. read only . init par block data */
10578 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
10579 /* ----------------------------------- */
10580 /* The binomial coefficients form a triangular matrix. */
10581 /* This matrix is completed in table CNP by its transposition. */
10582 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
10584 /* Initialisation is done by block-data MMLLL09.RES, */
10585 /* created by program MQINICNP.FOR (see the team (AC) ). */
10587 /* **********************************************************************
10592 /* ***********************************************************************
10595 /* Parameter adjustments */
10599 /* Function Body */
10600 if (*ncoeff < 1 || *ncoeff - 1 > 60) {
10606 /* CONSTANT TERM OF THE NEW CURVE */
10612 for (k = 2; k <= i__1; ++k) {
10613 cij1 += crvold[k * 3 + 1];
10614 cij2 += crvold[k * 3 + 2];
10615 cij3 += crvold[k * 3 + 3];
10621 if (*ncoeff == 1) {
10625 /* INTERMEDIARY POWER OF THE PARAMETER */
10627 ncfm1 = *ncoeff - 1;
10630 for (j = 2; j <= i__1; ++j) {
10632 cij1 = crvold[j * 3 + 1];
10633 cij2 = crvold[j * 3 + 2];
10634 cij3 = crvold[j * 3 + 3];
10636 for (k = j + 1; k <= i__2; ++k) {
10637 bid = mmcmcnp_.cnp[k - 1 + (j - 1) * 61];
10638 cij1 += crvold[k * 3 + 1] * bid;
10639 cij2 += crvold[k * 3 + 2] * bid;
10640 cij3 += crvold[k * 3 + 3] * bid;
10643 crvnew[j * 3 + 1] = cij1 * m1jm1;
10644 crvnew[j * 3 + 2] = cij2 * m1jm1;
10645 crvnew[j * 3 + 3] = cij3 * m1jm1;
10649 /* TERM OF THE HIGHEST DEGREE */
10651 crvnew[*ncoeff * 3 + 1] = -crvold[*ncoeff * 3 + 1] * m1jm1;
10652 crvnew[*ncoeff * 3 + 2] = -crvold[*ncoeff * 3 + 2] * m1jm1;
10653 crvnew[*ncoeff * 3 + 3] = -crvold[*ncoeff * 3 + 3] * m1jm1;
10656 AdvApp2Var_SysBase::maermsg_("MVCVINV", iercod, 7L);
10660 //=======================================================================
10661 //function : mvgaus0_
10663 //=======================================================================
10664 int mvgaus0_(integer *kindic,
10665 doublereal *urootl,
10666 doublereal *hiltab,
10671 /* System generated locals */
10674 /* Local variables */
10675 doublereal tamp[40];
10676 integer ndegl, kg, ii;
10678 /* **********************************************************************
10683 /* Loading of a degree gives roots of LEGENDRE polynom */
10684 /* DEFINED on [-1,1] and weights of Gauss quadrature formulas */
10685 /* (based on corresponding LAGRANGIAN interpolators). */
10686 /* The symmetry relative to 0 is used between [-1,0] and [0,1]. */
10690 /* . VOLUMIC, LEGENDRE, LAGRANGE, GAUSS */
10692 /* INPUT ARGUMENTSE : */
10693 /* ------------------ */
10695 /* KINDIC : Takes values from 1 to 10 depending of the degree */
10696 /* of the used polynom. */
10697 /* The degree of the polynom is equal to 4 k, i.e. 4, 8, */
10698 /* 12, 16, 20, 24, 28, 32, 36 and 40. */
10700 /* OUTPUT ARGUMENTS : */
10701 /* ------------------- */
10703 /* UROOTL : Roots of LEGENDRE polynom in domain [1,0] */
10704 /* given in decreasing order. For domain [-1,0], it is */
10705 /* necessary to take the opposite values. */
10706 /* HILTAB : LAGRANGE interpolators associated to roots. For */
10707 /* opposed roots, interpolatorsare equal. */
10708 /* NBRVAL : Nb of coefficients. Is equal to the half of degree */
10709 /* depending on the symmetry (i.e. 2*KINDIC). */
10711 /* IERCOD : Error code: */
10712 /* < 0 ==> Attention - Warning */
10713 /* =-1 ==> Value of false KINDIC. NBRVAL is forced to 20 */
10715 /* = 0 ==> Everything is OK */
10717 /* COMMON USED : */
10718 /* ---------------- */
10720 /* REFERENCES CALLED : */
10721 /* ------------------- */
10723 /* DESCRIPTION/NOTES/LIMITATIONS : */
10724 /* --------------------------------- */
10725 /* If KINDIC is not correct (i.e < 1 or > 10), the degree is set */
10726 /* to 40 directly (ATTENTION to overload - to avoid it, */
10727 /* preview UROOTL and HILTAB dimensioned at least to 20). */
10729 /* The value of coefficients was calculated with quadruple precision
10730 /* by JJM with help of GD. */
10731 /* Checking of roots was done by GD. */
10733 /* See detailed explications on the listing */
10735 /* **********************************************************************
10739 /* ------------------------------------ */
10740 /* ****** Test validity of KINDIC ** */
10741 /* ------------------------------------ */
10743 /* Parameter adjustments */
10747 /* Function Body */
10750 if (kg < 1 || kg > 10) {
10755 ndegl = *nbrval << 1;
10757 /* ----------------------------------------------------------------------
10759 /* ****** Load NBRVAL positive roots depending on the degree **
10761 /* ----------------------------------------------------------------------
10763 /* ATTENTION : Sign minus (-) in the loop is intentional. */
10765 mmextrl_(&ndegl, tamp);
10767 for (ii = 1; ii <= i__1; ++ii) {
10768 urootl[ii] = -tamp[ii - 1];
10772 /* ------------------------------------------------------------------- */
10773 /* ****** Loading of NBRVAL Gauss weight depending on the degree ** */
10774 /* ------------------------------------------------------------------- */
10776 mmexthi_(&ndegl, tamp);
10778 for (ii = 1; ii <= i__1; ++ii) {
10779 hiltab[ii] = tamp[ii - 1];
10783 /* ------------------------------- */
10784 /* ****** End of sub-program ** */
10785 /* ------------------------------- */
10790 //=======================================================================
10791 //function : mvpscr2_
10793 //=======================================================================
10794 int mvpscr2_(integer *ncoeff,
10795 doublereal *curve2,
10796 doublereal *tparam,
10797 doublereal *pntcrb)
10799 /* System generated locals */
10802 /* Local variables */
10804 doublereal xxx, yyy;
10808 /* **********************************************************************
10813 /* POSITIONING ON CURVE (NCF,2) IN SPACE OF DIMENSION 2. */
10817 /* TOUS,MATH_ACCES:: COURBE&,POSITIONNEMENT,&POINT. */
10819 /* INPUT ARGUMENTS : */
10820 /* ------------------ */
10821 /* NCOEFF : NUMBER OF COEFFICIENTS OF THE CURVE */
10822 /* CURVE2 : EQUATION OF CURVE 2D */
10823 /* TPARAM : VALUE OF PARAMETER AT GIVEN POINT */
10825 /* OUTPUT ARGUMENTS : */
10826 /* ------------------- */
10827 /* PNTCRB : COORDINATES OF POINT CORRESPONDING TO PARAMETER */
10828 /* TPARAM ON CURVE 2D CURVE2. */
10830 /* COMMONS USED : */
10831 /* ---------------- */
10833 /* REFERENCES CALLED : */
10834 /* ---------------------- */
10836 /* DESCRIPTION/NOTES/LIMITATIONS : */
10837 /* ----------------------------------- */
10838 /* MSCHEMA OF HORNER. */
10841 /* **********************************************************************
10845 /* -------- INITIALIZATIONS AND PROCESSING OF PARTICULAR CASES ----------
10848 /* ---> Cas when NCOEFF > 1 (case STANDARD). */
10849 /* Parameter adjustments */
10853 /* Function Body */
10854 if (*ncoeff >= 2) {
10857 /* ---> Case when NCOEFF <= 1. */
10858 if (*ncoeff <= 0) {
10862 } else if (*ncoeff == 1) {
10863 pntcrb[1] = curve2[3];
10864 pntcrb[2] = curve2[4];
10868 /* -------------------- MSCHEMA OF HORNER (PARTICULAR CASE) --------------
10873 if (*tparam == 1.) {
10877 for (kk = 1; kk <= i__1; ++kk) {
10878 xxx += curve2[(kk << 1) + 1];
10879 yyy += curve2[(kk << 1) + 2];
10883 } else if (*tparam == 0.) {
10884 pntcrb[1] = curve2[3];
10885 pntcrb[2] = curve2[4];
10889 /* ---------------------------- MSCHEMA OF HORNER ------------------------
10891 /* ---> TPARAM is different from 1.D0 and 0.D0. */
10893 ndeg = *ncoeff - 1;
10894 xxx = curve2[(*ncoeff << 1) + 1];
10895 yyy = curve2[(*ncoeff << 1) + 2];
10896 for (kk = ndeg; kk >= 1; --kk) {
10897 xxx = xxx * *tparam + curve2[(kk << 1) + 1];
10898 yyy = yyy * *tparam + curve2[(kk << 1) + 2];
10903 /* ------------------------ RECOVER THE CALCULATED POINT ---------------
10910 /* ------------------------------ THE END -------------------------------
10917 //=======================================================================
10918 //function : mvpscr3_
10920 //=======================================================================
10921 int mvpscr3_(integer *ncoeff,
10922 doublereal *curve3,
10923 doublereal *tparam,
10924 doublereal *pntcrb)
10927 /* System generated locals */
10930 /* Local variables */
10932 doublereal xxx, yyy, zzz;
10936 /* **********************************************************************
10941 /* POSITIONING ON A CURVE (3,NCF) IN THE SPACE OF DIMENSION 3. */
10945 /* TOUS, MATH_ACCES:: COURBE&,POSITIONNEMENT,&POINT. */
10947 /* INPUT ARGUMENTS : */
10948 /* ------------------ */
10949 /* NCOEFF : NB OF COEFFICIENTS OF THE CURVE */
10950 /* CURVE3 : EQUATION OF CURVE 3D */
10951 /* TPARAM : VALUE OF THE PARAMETER AT THE GIVEN POINT */
10953 /* OUTPUT ARGUMENTS : */
10954 /* ------------------- */
10955 /* PNTCRB : COORDINATES OF THE POINT CORRESPONDING TO PARAMETER */
10956 /* TPARAM ON CURVE 3D CURVE3. */
10958 /* COMMONS USED : */
10959 /* ---------------- */
10961 /* REFERENCES CALLED : */
10962 /* ---------------------- */
10965 /* DESCRIPTION/NOTES/LIMITATIONS : */
10966 /* ----------------------------------- */
10967 /* MSCHEMA OF HORNER. */
10969 /* **********************************************************************
10972 /* **********************************************************************
10976 /* -------- INITIALISATIONS AND PROCESSING OF PARTICULAR CASES ----------
10979 /* ---> Case when NCOEFF > 1 (cas STANDARD). */
10980 /* Parameter adjustments */
10984 /* Function Body */
10985 if (*ncoeff >= 2) {
10988 /* ---> Case when NCOEFF <= 1. */
10989 if (*ncoeff <= 0) {
10994 } else if (*ncoeff == 1) {
10995 pntcrb[1] = curve3[4];
10996 pntcrb[2] = curve3[5];
10997 pntcrb[3] = curve3[6];
11001 /* -------------------- MSCHEMA OF HORNER (PARTICULAR CASE) --------------
11006 if (*tparam == 1.) {
11011 for (kk = 1; kk <= i__1; ++kk) {
11012 xxx += curve3[kk * 3 + 1];
11013 yyy += curve3[kk * 3 + 2];
11014 zzz += curve3[kk * 3 + 3];
11018 } else if (*tparam == 0.) {
11019 pntcrb[1] = curve3[4];
11020 pntcrb[2] = curve3[5];
11021 pntcrb[3] = curve3[6];
11025 /* ---------------------------- MSCHEMA OF HORNER ------------------------
11027 /* ---> Here TPARAM is different from 1.D0 and 0.D0. */
11029 ndeg = *ncoeff - 1;
11030 xxx = curve3[*ncoeff * 3 + 1];
11031 yyy = curve3[*ncoeff * 3 + 2];
11032 zzz = curve3[*ncoeff * 3 + 3];
11033 for (kk = ndeg; kk >= 1; --kk) {
11034 xxx = xxx * *tparam + curve3[kk * 3 + 1];
11035 yyy = yyy * *tparam + curve3[kk * 3 + 2];
11036 zzz = zzz * *tparam + curve3[kk * 3 + 3];
11041 /* ------------------------ RETURN THE CALCULATED POINT ------------------
11049 /* ------------------------------ THE END -------------------------------
11056 //=======================================================================
11057 //function : AdvApp2Var_MathBase::mvsheld_
11059 //=======================================================================
11060 int AdvApp2Var_MathBase::mvsheld_(integer *n,
11066 /* System generated locals */
11067 integer dtab_dim1, dtab_offset, i__1, i__2;
11069 /* Local variables */
11072 integer i3, i4, i5, incrp1;
11075 /************************************************************************
11080 /* PARSING OF COLUMNS OF TABLE OF REAL*8 BY SHELL METHOD*/
11081 /* (IN INCREASING ORDER) */
11085 /* POINT-ENTRY, PARSING, SHELL */
11087 /* INPUT ARGUMENTS : */
11088 /* ------------------ */
11089 /* N : NUMBER OF COLUMNS OF THE TABLE */
11090 /* IS : NUMBER OF LINE OF THE TABLE */
11091 /* DTAB : TABLE OF REAL*8 TO BE PARSED */
11092 /* ICLE : POSITION OF THE KEY ON THE COLUMN */
11094 /* OUTPUT ARGUMENTS : */
11095 /* ------------------- */
11096 /* DTAB : PARSED TABLE */
11098 /* COMMONS USED : */
11099 /* ---------------- */
11102 /* REFERENCES CALLED : */
11103 /* ---------------------- */
11106 /* DESCRIPTION/NOTES/LIMITATIONS : */
11107 /* ----------------------------------- */
11108 /* CLASSIC SHELL METHOD : PARSING BY SERIES */
11109 /* Declaration DTAB(IS, 1) corresponds to DTAB(IS, *) */
11111 /* ***********************************************************************
11115 /* Parameter adjustments */
11117 dtab_offset = dtab_dim1 + 1;
11118 dtab -= dtab_offset;
11120 /* Function Body */
11124 /* ------------------------ */
11126 /* INITIALIZATION OF THE SEQUENCE OF INCREMENTS */
11127 /* FIND THE GREATEST INCREMENT SO THAT INCR < N/9 */
11131 if (incr >= *n / 9) {
11134 /* ----------------------------- */
11135 incr = incr * 3 + 1;
11138 /* LOOP ON INCREMENTS TILL INCR = 1 */
11139 /* PARSING BY SERIES DISTANT FROM INCR */
11143 /* ----------------- */
11145 for (i3 = incrp1; i3 <= i__1; ++i3) {
11146 /* ---------------------- */
11148 /* SET ELEMENT I3 AT ITS PLACE IN THE SERIES */
11155 /* ------------------------- */
11156 if (dtab[*icle + i4 * dtab_dim1] <= dtab[*icle + (i4 + incr) *
11162 for (i5 = 1; i5 <= i__2; ++i5) {
11163 /* ------------------ */
11164 dsave = dtab[i5 + i4 * dtab_dim1];
11165 dtab[i5 + i4 * dtab_dim1] = dtab[i5 + (i4 + incr) * dtab_dim1];
11166 dtab[i5 + (i4 + incr) * dtab_dim1] = dsave;
11177 /* PASSAGE TO THE NEXT INCREMENT */
11188 //=======================================================================
11189 //function : AdvApp2Var_MathBase::mzsnorm_
11191 //=======================================================================
11192 doublereal AdvApp2Var_MathBase::mzsnorm_(integer *ndimen,
11193 doublereal *vecteu)
11196 /* System generated locals */
11198 doublereal ret_val, d__1, d__2;
11200 /* Local variables */
11202 integer i__, irmax;
11206 /* ***********************************************************************
11211 /* SERVES to calculate the euclidian norm of a vector : */
11212 /* ____________________________ */
11213 /* Z = V V(1)**2 + V(2)**2 + ... */
11219 /* INPUT ARGUMENTS : */
11220 /* ------------------ */
11221 /* NDIMEN : Dimension of the vector */
11222 /* VECTEU : vector of dimension NDIMEN */
11224 /* OUTPUT ARGUMENTS : */
11225 /* ------------------- */
11226 /* MZSNORM : Value of the euclidian norm of vector VECTEU */
11228 /* COMMONS USED : */
11229 /* ---------------- */
11233 /* REFERENCES CALLED : */
11234 /* ---------------------- */
11236 /* R*8 ABS R*8 SQRT */
11238 /* DESCRIPTION/NOTESS/LIMITATIONS : */
11239 /* ----------------------------------- */
11240 /* To limit the risks of overflow, */
11241 /* the term of the strongest absolute value is factorized : */
11242 /* _______________________ */
11243 /* Z = !V(1)! * V 1 + (V(2)/V(1))**2 + ... */
11246 /* ***********************************************************************
11249 /* ***********************************************************************
11253 /* ***********************************************************************
11256 /* ***********************************************************************
11259 /* ___ Find the strongest absolute value term */
11261 /* Parameter adjustments */
11264 /* Function Body */
11267 for (i__ = 2; i__ <= i__1; ++i__) {
11268 if ((d__1 = vecteu[irmax], advapp_abs(d__1)) < (d__2 = vecteu[i__], advapp_abs(d__2)
11275 /* ___ Calculate the norme */
11277 if ((d__1 = vecteu[irmax], advapp_abs(d__1)) < 1.) {
11280 for (i__ = 1; i__ <= i__1; ++i__) {
11281 /* Computing 2nd power */
11282 d__1 = vecteu[i__];
11283 xsom += d__1 * d__1;
11286 ret_val = sqrt(xsom);
11290 for (i__ = 1; i__ <= i__1; ++i__) {
11291 if (i__ == irmax) {
11294 /* Computing 2nd power */
11295 d__1 = vecteu[i__] / vecteu[irmax];
11296 xsom += d__1 * d__1;
11300 ret_val = (d__1 = vecteu[irmax], advapp_abs(d__1)) * sqrt(xsom);
11303 /* ***********************************************************************
11305 /* RETURN CALLING PROGRAM */
11306 /* ***********************************************************************