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,
269 static integer c__8 = 8;
270 /* System generated locals */
274 /* Local variables */
276 static doublereal differ[100];
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 /* ***********************************************************************
335 AdvApp2Var_SysBase::mcrrqst_(&c__8, ndimen, differ, &iofset, &ier);
338 /* --- If allocation is refused, the trivial method is applied. */
344 for (i__ = 1; i__ <= i__1; ++i__) {
345 /* Computing 2nd power */
346 d__1 = point1[i__] - point2[i__];
347 *distan += d__1 * d__1;
349 *distan = sqrt(*distan);
351 /* --- Otherwise MZSNORM is used to minimize the risks of overflow
356 for (i__ = 1; i__ <= i__1; ++i__) {
358 differ[j] = point2[i__] - point1[i__];
361 *distan = AdvApp2Var_MathBase::mzsnorm_(ndimen, &differ[iofset]);
365 /* ***********************************************************************
367 /* RETURN CALLING PROGRAM */
368 /* ***********************************************************************
371 /* --- Dynamic Desallocation */
374 AdvApp2Var_SysBase::mcrdelt_(&c__8, ndimen, differ, &iofset, &ier);
380 //=======================================================================
383 //=======================================================================
384 int mfac_(doublereal *f,
388 /* System generated locals */
391 /* Local variables */
394 /* FORTRAN CONFORME AU TEXT */
395 /* CALCUL DE MFACTORIEL N */
396 /* Parameter adjustments */
402 for (i__ = 2; i__ <= i__1; ++i__) {
404 f[i__] = i__ * f[i__ - 1];
409 //=======================================================================
410 //function : AdvApp2Var_MathBase::mmapcmp_
412 //=======================================================================
413 int AdvApp2Var_MathBase::mmapcmp_(integer *ndim,
420 /* System generated locals */
421 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
424 /* Local variables */
425 static integer ipair, nd, ndegre, impair, ibb, idg;
426 //extern int mgsomsg_();//mgenmsg_(),
428 /* **********************************************************************
433 /* Compression of curve CRVOLD in a table of */
434 /* coeff. of even : CRVNEW(*,0,*) */
435 /* and uneven range : CRVNEW(*,1,*). */
439 /* COMPRESSION,CURVE. */
441 /* INPUT ARGUMENTS : */
442 /* ------------------ */
443 /* NDIM : Space Dimension. */
444 /* NCOFMX : Max nb of coeff. of the curve to compress. */
445 /* NCOEFF : Max nb of coeff. of the compressed curve. */
446 /* CRVOLD : The curve (0:NCOFMX-1,NDIM) to compress. */
448 /* OUTPUT ARGUMENTS : */
449 /* ------------------- */
450 /* CRVNEW : Curve compacted in (0:(NCOEFF-1)/2,0,NDIM) (containing
452 /* even terms) and in (0:(NCOEFF-1)/2,1,NDIM) */
453 /* (containing uneven terms). */
456 /* ---------------- */
458 /* REFERENCES CALLED : */
459 /* ----------------------- */
461 /* DESCRIPTION/NOTES/LIMITATIONS : */
462 /* ----------------------------------- */
463 /* This routine is useful to prepare coefficients of a */
464 /* curve in an orthogonal base (Legendre or Jacobi) before */
465 /* calculating the coefficients in the canonical; base [-1,1] by */
467 /* ***********************************************************************
470 /* Name of the routine */
472 /* Parameter adjustments */
473 crvold_dim1 = *ncofmx;
474 crvold_offset = crvold_dim1;
475 crvold -= crvold_offset;
476 crvnew_dim1 = (*ncoeff - 1) / 2 + 1;
477 crvnew_offset = crvnew_dim1 << 1;
478 crvnew -= crvnew_offset;
481 ibb = AdvApp2Var_SysBase::mnfndeb_();
483 AdvApp2Var_SysBase::mgenmsg_("MMAPCMP", 7L);
486 ndegre = *ncoeff - 1;
488 for (nd = 1; nd <= i__1; ++nd) {
491 for (idg = 0; idg <= i__2; ++idg) {
492 crvnew[idg + (nd << 1) * crvnew_dim1] = crvold[ipair + nd *
501 i__2 = (ndegre - 1) / 2;
502 for (idg = 0; idg <= i__2; ++idg) {
503 crvnew[idg + ((nd << 1) + 1) * crvnew_dim1] = crvold[impair + nd *
514 /* ---------------------------------- The end ---------------------------
518 AdvApp2Var_SysBase::mgsomsg_("MMAPCMP", 7L);
523 //=======================================================================
524 //function : mmaper0_
526 //=======================================================================
527 int mmaper0_(integer *ncofmx,
536 /* System generated locals */
537 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
540 /* Local variables */
542 static doublereal bidon;
543 static integer ii, nd;
545 /* ***********************************************************************
550 /* Calculate the max error of approximation done when */
551 /* only the first NCFNEW coefficients of a curve are preserved.
553 /* Degree NCOEFF-1 written in the base of Legendre (Jacobi */
558 /* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
560 /* INPUT ARGUMENTS : */
561 /* ------------------ */
562 /* NCOFMX : Max. degree of the curve. */
563 /* NDIMEN : Space dimension. */
564 /* NCOEFF : Degree +1 of the curve. */
565 /* CRVLGD : Curve the degree which of should be lowered. */
566 /* NCFNEW : Degree +1 of the resulting polynom. */
568 /* OUTPUT ARGUMENTS : */
569 /* ------------------- */
570 /* YCVMAX : Auxiliary Table (max error on each dimension).
572 /* ERRMAX : Precision of the approximation. */
575 /* ---------------- */
577 /* REFERENCES CALLED : */
578 /* ----------------------- */
580 /* DESCRIPTION/NOTES/LIMITATIONS : */
581 /* ----------------------------------- */
582 /* ***********************************************************************
586 /* ------------------- Init to calculate an error -----------------------
589 /* Parameter adjustments */
591 crvlgd_dim1 = *ncofmx;
592 crvlgd_offset = crvlgd_dim1 + 1;
593 crvlgd -= crvlgd_offset;
597 for (ii = 1; ii <= i__1; ++ii) {
602 /* ------ Minimum that can be reached : Stop at 1 or NCFNEW ------
606 if (*ncfnew + 1 > ncut) {
610 /* -------------- Elimination of high degree coefficients-----------
612 /* ----------- Loop on the series of Legendre: NCUT --> NCOEFF --------
616 for (ii = ncut; ii <= i__1; ++ii) {
617 /* Factor of renormalization (Maximum of Li(t)). */
618 bidon = ((ii - 1) * 2. + 1.) / 2.;
622 for (nd = 1; nd <= i__2; ++nd) {
623 ycvmax[nd] += (d__1 = crvlgd[ii + nd * crvlgd_dim1], advapp_abs(d__1)) *
630 /* -------------- The error is the norm of the vector error ---------------
633 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
635 /* --------------------------------- Fin --------------------------------
641 //=======================================================================
642 //function : mmaper2_
644 //=======================================================================
645 int mmaper2_(integer *ncofmx,
654 /* Initialized data */
656 static doublereal xmaxj[57] = { .9682458365518542212948163499456,
657 .986013297183269340427888048593603,
658 1.07810420343739860362585159028115,
659 1.17325804490920057010925920756025,
660 1.26476561266905634732910520370741,
661 1.35169950227289626684434056681946,
662 1.43424378958284137759129885012494,
663 1.51281316274895465689402798226634,
664 1.5878364329591908800533936587012,
665 1.65970112228228167018443636171226,
666 1.72874345388622461848433443013543,
667 1.7952515611463877544077632304216,
668 1.85947199025328260370244491818047,
669 1.92161634324190018916351663207101,
670 1.98186713586472025397859895825157,
671 2.04038269834980146276967984252188,
672 2.09730119173852573441223706382076,
673 2.15274387655763462685970799663412,
674 2.20681777186342079455059961912859,
675 2.25961782459354604684402726624239,
676 2.31122868752403808176824020121524,
677 2.36172618435386566570998793688131,
678 2.41117852396114589446497298177554,
679 2.45964731268663657873849811095449,
680 2.50718840313973523778244737914028,
681 2.55385260994795361951813645784034,
682 2.59968631659221867834697883938297,
683 2.64473199258285846332860663371298,
684 2.68902863641518586789566216064557,
685 2.73261215675199397407027673053895,
686 2.77551570192374483822124304745691,
687 2.8177699459714315371037628127545,
688 2.85940333797200948896046563785957,
689 2.90044232019793636101516293333324,
690 2.94091151970640874812265419871976,
691 2.98083391718088702956696303389061,
692 3.02023099621926980436221568258656,
693 3.05912287574998661724731962377847,
694 3.09752842783622025614245706196447,
695 3.13546538278134559341444834866301,
696 3.17295042316122606504398054547289,
697 3.2099992681699613513775259670214,
698 3.24662674946606137764916854570219,
699 3.28284687953866689817670991319787,
700 3.31867291347259485044591136879087,
701 3.35411740487202127264475726990106,
702 3.38919225660177218727305224515862,
703 3.42390876691942143189170489271753,
704 3.45827767149820230182596660024454,
705 3.49230918177808483937957161007792,
706 3.5260130200285724149540352829756,
707 3.55939845146044235497103883695448,
708 3.59247431368364585025958062194665,
709 3.62524904377393592090180712976368,
710 3.65773070318071087226169680450936,
711 3.68992700068237648299565823810245,
712 3.72184531357268220291630708234186 };
714 /* System generated locals */
715 integer crvjac_dim1, crvjac_offset, i__1, i__2;
718 /* Local variables */
719 static integer idec, ncut;
720 static doublereal bidon;
721 static integer ii, nd;
725 /* ***********************************************************************
730 /* Calculate max approximation error i faite lorsque l' on */
731 /* ne conserve que les premiers NCFNEW coefficients d' une courbe
733 /* de degre NCOEFF-1 ecrite dans la base de Jacobi d' ordre 2. */
737 /* JACOBI, POLYGON, APPROXIMATION, ERROR. */
739 /* INPUT ARGUMENTS : */
740 /* ------------------ */
741 /* NCOFMX : Max. degree of the curve. */
742 /* NDIMEN : Space dimension. */
743 /* NCOEFF : Degree +1 of the curve. */
744 /* CRVLGD : Curve the degree which of should be lowered. */
745 /* NCFNEW : Degree +1 of the resulting polynom. */
747 /* OUTPUT ARGUMENTS : */
748 /* ------------------- */
749 /* YCVMAX : Auxiliary Table (max error on each dimension).
751 /* ERRMAX : Precision of the approximation. */
754 /* ---------------- */
756 /* REFERENCES CALLED : */
757 /* ----------------------- */
758 /* DESCRIPTION/NOTES/LIMITATIONS : */
759 /* ----------------------------------- */
763 /* ------------------ Table of maximums of (1-t2)*Ji(t) ----------------
766 /* Parameter adjustments */
768 crvjac_dim1 = *ncofmx;
769 crvjac_offset = crvjac_dim1 + 1;
770 crvjac -= crvjac_offset;
776 /* ------------------- Init for error calculation -----------------------
780 for (ii = 1; ii <= i__1; ++ii) {
785 /* ------ Min. Degree that can be attained : Stop at 3 or NCFNEW ------
790 i__1 = idec, i__2 = *ncfnew + 1;
791 ncut = advapp_max(i__1,i__2);
793 /* -------------- Removal of coefficients of high degree -----------
795 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
799 for (ii = ncut; ii <= i__1; ++ii) {
800 /* Factor of renormalization. */
801 bidon = xmaxj[ii - idec];
803 for (nd = 1; nd <= i__2; ++nd) {
804 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
811 /* -------------- The error is the norm of the vector error ---------------
814 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
816 /* --------------------------------- Fin --------------------------------
822 /* MAPER4.f -- translated by f2c (version 19960827).
823 You must link the resulting object file with the libraries:
824 -lf2c -lm (in that order)
828 //=======================================================================
829 //function : mmaper4_
831 //=======================================================================
832 int mmaper4_(integer *ncofmx,
840 /* Initialized data */
842 static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
843 1.05299572648705464724876659688996,
844 1.0949715351434178709281698645813,
845 1.15078388379719068145021100764647,
846 1.2094863084718701596278219811869,
847 1.26806623151369531323304177532868,
848 1.32549784426476978866302826176202,
849 1.38142537365039019558329304432581,
850 1.43575531950773585146867625840552,
851 1.48850442653629641402403231015299,
852 1.53973611681876234549146350844736,
853 1.58953193485272191557448229046492,
854 1.63797820416306624705258190017418,
855 1.68515974143594899185621942934906,
856 1.73115699602477936547107755854868,
857 1.77604489805513552087086912113251,
858 1.81989256661534438347398400420601,
859 1.86276344480103110090865609776681,
860 1.90471563564740808542244678597105,
861 1.94580231994751044968731427898046,
862 1.98607219357764450634552790950067,
863 2.02556989246317857340333585562678,
864 2.06433638992049685189059517340452,
865 2.10240936014742726236706004607473,
866 2.13982350649113222745523925190532,
867 2.17661085564771614285379929798896,
868 2.21280102016879766322589373557048,
869 2.2484214321456956597803794333791,
870 2.28349755104077956674135810027654,
871 2.31805304852593774867640120860446,
872 2.35210997297725685169643559615022,
873 2.38568889602346315560143377261814,
874 2.41880904328694215730192284109322,
875 2.45148841120796359750021227795539,
876 2.48374387161372199992570528025315,
877 2.5155912654873773953959098501893,
878 2.54704548720896557684101746505398,
879 2.57812056037881628390134077704127,
880 2.60882970619319538196517982945269,
881 2.63918540521920497868347679257107,
882 2.66919945330942891495458446613851,
883 2.69888301230439621709803756505788,
884 2.72824665609081486737132853370048,
885 2.75730041251405791603760003778285,
886 2.78605380158311346185098508516203,
887 2.81451587035387403267676338931454,
888 2.84269522483114290814009184272637,
889 2.87060005919012917988363332454033,
890 2.89823818258367657739520912946934,
891 2.92561704377132528239806135133273,
892 2.95274375377994262301217318010209,
893 2.97962510678256471794289060402033,
894 3.00626759936182712291041810228171,
895 3.03267744830655121818899164295959,
896 3.05886060707437081434964933864149 };
898 /* System generated locals */
899 integer crvjac_dim1, crvjac_offset, i__1, i__2;
902 /* Local variables */
903 static integer idec, ncut;
904 static doublereal bidon;
905 static integer ii, nd;
909 /* ***********************************************************************
914 /* Calculate the max. error of approximation made when */
915 /* only first NCFNEW coefficients of a curve are preserved
917 /* degree NCOEFF-1 is written in the base of Jacobi of order 4. */
920 /* LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
922 /* INPUT ARGUMENTS : */
923 /* ------------------ */
924 /* NCOFMX : Max. degree of the curve. */
925 /* NDIMEN : Space dimension. */
926 /* NCOEFF : Degree +1 of the curve. */
927 /* CRVJAC : Curve the degree which of should be lowered. */
928 /* NCFNEW : Degree +1 of the resulting polynom. */
930 /* OUTPUT ARGUMENTS : */
931 /* ------------------- */
932 /* YCVMAX : Auxiliary Table (max error on each dimension).
934 /* ERRMAX : Precision of the approximation. */
937 /* ---------------- */
939 /* REFERENCES CALLED : */
940 /* ----------------------- */
942 /* DESCRIPTION/NOTES/LIMITATIONS : */
945 /* ***********************************************************************
949 /* ---------------- Table of maximums of ((1-t2)2)*Ji(t) ---------------
952 /* Parameter adjustments */
954 crvjac_dim1 = *ncofmx;
955 crvjac_offset = crvjac_dim1 + 1;
956 crvjac -= crvjac_offset;
962 /* ------------------- Init for error calculation -----------------------
966 for (ii = 1; ii <= i__1; ++ii) {
971 /* ------ Min. Degree that can be attained : Stop at 5 or NCFNEW ------
976 i__1 = idec, i__2 = *ncfnew + 1;
977 ncut = advapp_max(i__1,i__2);
979 /* -------------- Removal of high degree coefficients -----------
981 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
985 for (ii = ncut; ii <= i__1; ++ii) {
986 /* Factor of renormalisation. */
987 bidon = xmaxj[ii - idec];
989 for (nd = 1; nd <= i__2; ++nd) {
990 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
997 /* -------------- The error is the norm of the error vector ---------------
1000 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
1002 /* --------------------------------- End --------------------------------
1008 //=======================================================================
1009 //function : mmaper6_
1011 //=======================================================================
1012 int mmaper6_(integer *ncofmx,
1021 /* Initialized data */
1023 static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
1024 1.11626917091567929907256116528817,
1025 1.1327140810290884106278510474203,
1026 1.1679452722668028753522098022171,
1027 1.20910611986279066645602153641334,
1028 1.25228283758701572089625983127043,
1029 1.29591971597287895911380446311508,
1030 1.3393138157481884258308028584917,
1031 1.3821288728999671920677617491385,
1032 1.42420414683357356104823573391816,
1033 1.46546895108549501306970087318319,
1034 1.50590085198398789708599726315869,
1035 1.54550385142820987194251585145013,
1036 1.58429644271680300005206185490937,
1037 1.62230484071440103826322971668038,
1038 1.65955905239130512405565733793667,
1039 1.69609056468292429853775667485212,
1040 1.73193098017228915881592458573809,
1041 1.7671112206990325429863426635397,
1042 1.80166107681586964987277458875667,
1043 1.83560897003644959204940535551721,
1044 1.86898184653271388435058371983316,
1045 1.90180515174518670797686768515502,
1046 1.93410285411785808749237200054739,
1047 1.96589749778987993293150856865539,
1048 1.99721027139062501070081653790635,
1049 2.02806108474738744005306947877164,
1050 2.05846864831762572089033752595401,
1051 2.08845055210580131460156962214748,
1052 2.11802334209486194329576724042253,
1053 2.14720259305166593214642386780469,
1054 2.17600297710595096918495785742803,
1055 2.20443832785205516555772788192013,
1056 2.2325216999457379530416998244706,
1057 2.2602654243075083168599953074345,
1058 2.28768115912702794202525264301585,
1059 2.3147799369092684021274946755348,
1060 2.34157220782483457076721300512406,
1061 2.36806787963276257263034969490066,
1062 2.39427635443992520016789041085844,
1063 2.42020656255081863955040620243062,
1064 2.44586699364757383088888037359254,
1065 2.47126572552427660024678584642791,
1066 2.49641045058324178349347438430311,
1067 2.52130850028451113942299097584818,
1068 2.54596686772399937214920135190177,
1069 2.5703922285006754089328998222275,
1070 2.59459096001908861492582631591134,
1071 2.61856915936049852435394597597773,
1072 2.64233265984385295286445444361827,
1073 2.66588704638685848486056711408168,
1074 2.68923766976735295746679957665724,
1075 2.71238965987606292679677228666411 };
1077 /* System generated locals */
1078 integer crvjac_dim1, crvjac_offset, i__1, i__2;
1081 /* Local variables */
1082 static integer idec, ncut;
1083 static doublereal bidon;
1084 static integer ii, nd;
1088 /* ***********************************************************************
1092 /* Calculate the max. error of approximation made when */
1093 /* only first NCFNEW coefficients of a curve are preserved
1095 /* degree NCOEFF-1 is written in the base of Jacobi of order 6. */
1098 /* JACOBI,POLYGON,APPROXIMATION,ERROR. */
1100 /* INPUT ARGUMENTS : */
1101 /* ------------------ */
1102 /* NCOFMX : Max. degree of the curve. */
1103 /* NDIMEN : Space dimension. */
1104 /* NCOEFF : Degree +1 of the curve. */
1105 /* CRVJAC : Curve the degree which of should be lowered. */
1106 /* NCFNEW : Degree +1 of the resulting polynom. */
1108 /* OUTPUT ARGUMENTS : */
1109 /* ------------------- */
1110 /* YCVMAX : Auxiliary Table (max error on each dimension).
1112 /* ERRMAX : Precision of the approximation. */
1114 /* COMMONS USED : */
1115 /* ---------------- */
1117 /* REFERENCES CALLED : */
1118 /* ----------------------- */
1120 /* DESCRIPTION/NOTES/LIMITATIONS : */
1122 /* ***********************************************************************
1126 /* ---------------- Table of maximums of ((1-t2)3)*Ji(t) ---------------
1129 /* Parameter adjustments */
1131 crvjac_dim1 = *ncofmx;
1132 crvjac_offset = crvjac_dim1 + 1;
1133 crvjac -= crvjac_offset;
1139 /* ------------------- Init for error calculation -----------------------
1143 for (ii = 1; ii <= i__1; ++ii) {
1148 /* ------ Min Degree that can be attained : Stop at 3 or NCFNEW ------
1153 i__1 = idec, i__2 = *ncfnew + 1;
1154 ncut = advapp_max(i__1,i__2);
1156 /* -------------- Removal of high degree coefficients -----------
1158 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ----------
1162 for (ii = ncut; ii <= i__1; ++ii) {
1163 /* Factor of renormalization. */
1164 bidon = xmaxj[ii - idec];
1166 for (nd = 1; nd <= i__2; ++nd) {
1167 ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) *
1174 /* -------------- The error is the norm of the vector error ---------------
1177 *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
1179 /* --------------------------------- END --------------------------------
1185 //=======================================================================
1186 //function : AdvApp2Var_MathBase::mmaperx_
1188 //=======================================================================
1189 int AdvApp2Var_MathBase::mmaperx_(integer *ncofmx,
1200 /* System generated locals */
1201 integer crvjac_dim1, crvjac_offset;
1203 /* Local variables */
1204 static integer jord;
1206 /* **********************************************************************
1210 /* Calculate the max. error of approximation made when */
1211 /* only first NCFNEW coefficients of a curve are preserved
1213 /* degree NCOEFF-1 is written in the base of Jacobi of order IORDRE. */
1216 /* JACOBI,LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
1218 /* INPUT ARGUMENTS : */
1219 /* ------------------ */
1220 /* NCOFMX : Max. degree of the curve. */
1221 /* NDIMEN : Space dimension. */
1222 /* NCOEFF : Degree +1 of the curve. */
1223 /* IORDRE : Order of continuity at the extremities. */
1224 /* CRVJAC : Curve the degree which of should be lowered. */
1225 /* NCFNEW : Degree +1 of the resulting polynom. */
1227 /* OUTPUT ARGUMENTS : */
1228 /* ------------------- */
1229 /* YCVMAX : Auxiliary Table (max error on each dimension).
1231 /* ERRMAX : Precision of the approximation. */
1232 /* IERCOD = 0, OK */
1233 /* = 1, order of constraints (IORDRE) is not within the */
1234 /* autorized values. */
1235 /* COMMONS USED : */
1236 /* ---------------- */
1238 /* REFERENCES CALLED : */
1239 /* ----------------------- */
1241 /* DESCRIPTION/NOTES/LIMITATIONS : */
1242 /* ----------------------------------- */
1243 /* Canceled and replaced MMAPERR. */
1244 /* ***********************************************************************
1248 /* Parameter adjustments */
1250 crvjac_dim1 = *ncofmx;
1251 crvjac_offset = crvjac_dim1 + 1;
1252 crvjac -= crvjac_offset;
1256 /* --> Order of Jacobi polynoms */
1257 jord = ( *iordre + 1) << 1;
1260 mmaper0_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1262 } else if (jord == 2) {
1263 mmaper2_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1265 } else if (jord == 4) {
1266 mmaper4_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1268 } else if (jord == 6) {
1269 mmaper6_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1275 /* ----------------------------------- Fin ------------------------------
1281 //=======================================================================
1282 //function : mmarc41_
1284 //=======================================================================
1285 int mmarc41_(integer *ndimax,
1295 /* System generated locals */
1296 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
1299 /* Local variables */
1300 static integer nboct;
1301 static doublereal tbaux[61];
1303 static doublereal bid;
1304 static integer ncf, ncj;
1307 /* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
1308 /* IMPLICIT INTEGER (I-N) */
1310 /* ***********************************************************************
1315 /* Creation of curve C2(v) defined on (0,1) identic to */
1316 /* curve C1(u) defined on (U0,U1) (change of parameter */
1321 /* LIMITATION, RESTRICTION, CURVE */
1323 /* INPUT ARGUMENTS : */
1324 /* ------------------ */
1325 /* NDIMAX : Space Dimensioning. */
1326 /* NDIMEN : Curve Dimension. */
1327 /* NCOEFF : Nb of coefficients of the curve. */
1328 /* CRVOLD : Curve to be limited. */
1329 /* UPARA0 : Min limit of the interval limiting the curve.
1331 /* UPARA1 : Max limit of the interval limiting the curve.
1334 /* OUTPUT ARGUMENTS : */
1335 /* ------------------- */
1336 /* CRVNEW : Relimited curve, defined on (0,1) and equal to */
1337 /* CRVOLD defined on (U0,U1). */
1338 /* IERCOD : = 0, OK */
1339 /* =10, Nb of coeff. <1 or > 61. */
1341 /* COMMONS USED : */
1342 /* ---------------- */
1343 /* REFERENCES CALLED : */
1344 /* ---------------------- */
1346 /* MAERMSG MCRFILL MVCVIN2 */
1349 /* DESCRIPTION/NOTES/LIMITATIONS : */
1350 /* ----------------------------------- */
1351 /* ---> Algorithm used in this general case is based on the */
1352 /* following principle : */
1353 /* Let S(t) = a0 + a1*t + a2*t**2 + ... of degree NCOEFF-1, and */
1354 /* U(t) = b0 + b1*t, then the coeff. of */
1355 /* S(U(t)) are calculated step by step with help of table TBAUX. */
1356 /* At each step number N (N=2 to NCOEFF), TBAUX(n) contains */
1357 /* the n-th coefficient of U(t)**N for n=1 to N. (RBD) */
1358 /* ---> Reference : KNUTH, 'The Art of Computer Programming', */
1359 /* Vol. 2/'Seminumerical Algorithms', */
1360 /* Ex. 11 p:451 et solution p:562. (RBD) */
1362 /* ---> Removal of the input argument CRVOLD by CRVNEW is */
1363 /* possible, which means that the call : */
1364 /* CALL MMARC41(NDIMAX,NDIMEN,NCOEFF,CURVE,UPARA0,UPARA1 */
1365 /* ,CURVE,IERCOD) */
1366 /* is absolutely LEGAL. (RBD) */
1369 /* **********************************************************************
1372 /* Name of the routine */
1374 /* Auxiliary table of coefficients of (UPARA1-UPARA0)T+UPARA0 */
1375 /* with power N=1 to NCOEFF-1. */
1378 /* Parameter adjustments */
1379 crvnew_dim1 = *ndimax;
1380 crvnew_offset = crvnew_dim1 + 1;
1381 crvnew -= crvnew_offset;
1382 crvold_dim1 = *ndimax;
1383 crvold_offset = crvold_dim1 + 1;
1384 crvold -= crvold_offset;
1388 /* **********************************************************************
1390 /* CASE WHEN PROCESSING CAN'T BE DONE */
1391 /* **********************************************************************
1393 if (*ncoeff > 61 || *ncoeff < 1) {
1397 /* **********************************************************************
1400 /* **********************************************************************
1402 if (*ndimen == *ndimax && *upara0 == 0. && *upara1 == 1.) {
1403 nboct = (*ndimax << 3) * *ncoeff;
1404 AdvApp2Var_SysBase::mcrfill_((integer *)&nboct,
1405 (char *)&crvold[crvold_offset],
1406 (char *)&crvnew[crvnew_offset]);
1409 /* **********************************************************************
1411 /* INVERSION 3D : FAST PROCESSING */
1412 /* **********************************************************************
1414 if (*upara0 == 1. && *upara1 == 0.) {
1415 if (*ndimen == 3 && *ndimax == 3 && *ncoeff <= 21) {
1416 mvcvinv_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
1420 /* ******************************************************************
1422 /* INVERSION 2D : FAST PROCESSING */
1423 /* ******************************************************************
1425 if (*ndimen == 2 && *ndimax == 2 && *ncoeff <= 21) {
1426 mvcvin2_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset],
1431 /* **********************************************************************
1433 /* GENERAL PROCESSING */
1434 /* **********************************************************************
1436 /* -------------------------- Initializations ---------------------------
1440 for (nd = 1; nd <= i__1; ++nd) {
1441 crvnew[nd + crvnew_dim1] = crvold[nd + crvold_dim1];
1448 tbaux[1] = *upara1 - *upara0;
1450 /* ----------------------- Calculation of coeff. of CRVNEW ------------------
1454 for (ncf = 2; ncf <= i__1; ++ncf) {
1456 /* ------------ Take into account NCF-th coeff. of CRVOLD --------
1460 for (ncj = 1; ncj <= i__2; ++ncj) {
1461 bid = tbaux[ncj - 1];
1463 for (nd = 1; nd <= i__3; ++nd) {
1464 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + ncf *
1471 bid = tbaux[ncf - 1];
1473 for (nd = 1; nd <= i__2; ++nd) {
1474 crvnew[nd + ncf * crvnew_dim1] = crvold[nd + ncf * crvold_dim1] *
1479 /* --------- Calculate (NCF+1) coeff. of ((U1-U0)*t + U0)**(NCF) ---
1482 bid = *upara1 - *upara0;
1483 tbaux[ncf] = tbaux[ncf - 1] * bid;
1484 for (ncj = ncf; ncj >= 2; --ncj) {
1485 tbaux[ncj - 1] = tbaux[ncj - 1] * *upara0 + tbaux[ncj - 2] * bid;
1488 tbaux[0] *= *upara0;
1493 /* -------------- Take into account the last coeff. of CRVOLD -----------
1497 for (ncj = 1; ncj <= i__1; ++ncj) {
1498 bid = tbaux[ncj - 1];
1500 for (nd = 1; nd <= i__2; ++nd) {
1501 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + *ncoeff *
1508 for (nd = 1; nd <= i__1; ++nd) {
1509 crvnew[nd + *ncoeff * crvnew_dim1] = crvold[nd + *ncoeff *
1510 crvold_dim1] * tbaux[*ncoeff - 1];
1514 /* ---------------------------- The end ---------------------------------
1519 AdvApp2Var_SysBase::maermsg_("MMARC41", iercod, 7L);
1525 //=======================================================================
1526 //function : AdvApp2Var_MathBase::mmarcin_
1528 //=======================================================================
1529 int AdvApp2Var_MathBase::mmarcin_(integer *ndimax,
1539 /* System generated locals */
1540 integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1,
1544 /* Local variables */
1545 static doublereal x0, x1;
1547 static doublereal tabaux[61];
1549 static doublereal bid;
1552 static doublereal eps3;
1556 /* **********************************************************************
1559 /* Creation of curve C2(v) defined on [U0,U1] identic to */
1560 /* curve C1(u) defined on [-1,1] (change of parameter */
1561 /* of a curve) with INVERSION of indices of the resulting table. */
1565 /* GENERALIZED LIMITATION, RESTRICTION, INVERSION, CURVE */
1567 /* INPUT ARGUMENTS : */
1568 /* ------------------ */
1569 /* NDIMAX : Maximum Space Dimensioning. */
1570 /* NDIMEN : Curve Dimension. */
1571 /* NCOEFF : Nb of coefficients of the curve. */
1572 /* CRVOLD : Curve to be limited. */
1573 /* U0 : Min limit of the interval limiting the curve.
1575 /* U1 : Max limit of the interval limiting the curve.
1578 /* OUTPUT ARGUMENTS : */
1579 /* ------------------- */
1580 /* CRVNEW : Relimited curve, defined on [U0,U1] and equal to */
1581 /* CRVOLD defined on [-1,1]. */
1582 /* IERCOD : = 0, OK */
1583 /* =10, Nb of coeff. <1 or > 61. */
1584 /* =13, the requested interval of variation is null. */
1585 /* COMMONS USED : */
1586 /* ---------------- */
1587 /* REFERENCES CALLED : */
1588 /* ---------------------- */
1589 /* DESCRIPTION/NOTES/LIMITATIONS : */
1590 /* ----------------------------------- */
1592 /* **********************************************************************
1595 /* Name of the routine */
1597 /* Auxiliary table of coefficients of X1*T+X0 */
1598 /* with power N=1 to NCOEFF-1. */
1601 /* Parameter adjustments */
1602 crvnew_dim1 = *ndimax;
1603 crvnew_offset = crvnew_dim1 + 1;
1604 crvnew -= crvnew_offset;
1605 crvold_dim1 = *ncoeff;
1606 crvold_offset = crvold_dim1 + 1;
1607 crvold -= crvold_offset;
1610 ibb = AdvApp2Var_SysBase::mnfndeb_();
1612 AdvApp2Var_SysBase::mgenmsg_("MMARCIN", 7L);
1615 /* At zero machine it is tested if the output interval is not null */
1617 AdvApp2Var_MathBase::mmveps3_(&eps3);
1618 if ((d__1 = *u1 - *u0, advapp_abs(d__1)) < eps3) {
1624 /* **********************************************************************
1626 /* CASE WHEN THE PROCESSING IS IMPOSSIBLE */
1627 /* **********************************************************************
1629 if (*ncoeff > 61 || *ncoeff < 1) {
1633 /* **********************************************************************
1635 /* IF NO CHANGE OF THE INTERVAL OF DEFINITION */
1636 /* (ONLY INVERSION OF INDICES OF TABLE CRVOLD) */
1637 /* **********************************************************************
1639 if (*ndim == *ndimax && *u0 == -1. && *u1 == 1.) {
1640 AdvApp2Var_MathBase::mmcvinv_(ndim, ncoeff, ndim, &crvold[crvold_offset], &crvnew[
1644 /* **********************************************************************
1646 /* CASE WHEN THE NEW INTERVAL OF DEFINITION IS [0,1] */
1647 /* **********************************************************************
1649 if (*u0 == 0. && *u1 == 1.) {
1650 mmcvstd_(ncoeff, ndimax, ncoeff, ndim, &crvold[crvold_offset], &
1651 crvnew[crvnew_offset]);
1654 /* **********************************************************************
1656 /* GENERAL PROCESSING */
1657 /* **********************************************************************
1659 /* -------------------------- Initialization ---------------------------
1662 x0 = -(*u1 + *u0) / (*u1 - *u0);
1663 x1 = 2. / (*u1 - *u0);
1665 for (nd = 1; nd <= i__1; ++nd) {
1666 crvnew[nd + crvnew_dim1] = crvold[nd * crvold_dim1 + 1];
1675 /* ----------------------- Calculation of coeff. of CRVNEW ------------------
1679 for (ncf = 2; ncf <= i__1; ++ncf) {
1681 /* ------------ Take into account the NCF-th coeff. of CRVOLD --------
1685 for (ncj = 1; ncj <= i__2; ++ncj) {
1686 bid = tabaux[ncj - 1];
1688 for (nd = 1; nd <= i__3; ++nd) {
1689 crvnew[nd + ncj * crvnew_dim1] += crvold[ncf + nd *
1696 bid = tabaux[ncf - 1];
1698 for (nd = 1; nd <= i__2; ++nd) {
1699 crvnew[nd + ncf * crvnew_dim1] = crvold[ncf + nd * crvold_dim1] *
1704 /* --------- Calculation of (NCF+1) coeff. of [X1*t + X0]**(NCF) --------
1707 tabaux[ncf] = tabaux[ncf - 1] * x1;
1708 for (ncj = ncf; ncj >= 2; --ncj) {
1709 tabaux[ncj - 1] = tabaux[ncj - 1] * x0 + tabaux[ncj - 2] * x1;
1717 /* -------------- Take into account the last coeff. of CRVOLD -----------
1721 for (ncj = 1; ncj <= i__1; ++ncj) {
1722 bid = tabaux[ncj - 1];
1724 for (nd = 1; nd <= i__2; ++nd) {
1725 crvnew[nd + ncj * crvnew_dim1] += crvold[*ncoeff + nd *
1732 for (nd = 1; nd <= i__1; ++nd) {
1733 crvnew[nd + *ncoeff * crvnew_dim1] = crvold[*ncoeff + nd *
1734 crvold_dim1] * tabaux[*ncoeff - 1];
1738 /* ---------------------------- The end ---------------------------------
1743 AdvApp2Var_SysBase::maermsg_("MMARCIN", iercod, 7L);
1746 AdvApp2Var_SysBase::mgsomsg_("MMARCIN", 7L);
1751 //=======================================================================
1752 //function : mmatvec_
1754 //=======================================================================
1755 int mmatvec_(integer *nligne,
1766 /* System generated locals */
1769 /* Local variables */
1770 static logical ldbg;
1771 static integer jmin, jmax, i__, j, k;
1772 static doublereal somme;
1776 /* ***********************************************************************
1781 /* Produce vector matrix in form of profile */
1786 /* RESERVE, MATRIX, PRODUCT, VECTOR, PROFILE */
1788 /* INPUT ARGUMENTS : */
1789 /* -------------------- */
1790 /* NLIGNE : Line number of the matrix of constraints */
1791 /* NCOLON : Number of column of the matrix of constraints */
1792 /* GNSTOC: Number of coefficients in the profile of matrix GMATRI */
1794 /* GPOSIT: Table of positioning of terms of storage */
1795 /* GPOSIT(1,I) contains the number of terms-1 on the line I
1796 /* in the profile of the matrix. */
1797 /* GPOSIT(2,I) contains the index of storage of diagonal term*/
1799 /* GPOSIT(3,I) contains the index of column of the first term of */
1800 /* profile of line I */
1801 /* GNSTOC: Number of coefficients in the profile of matrix */
1803 /* GMATRI : Matrix of constraints in form of profile */
1804 /* VECIN : Input vector */
1805 /* DEBLIG : Line indexusing which the vector matrix is calculated */
1807 /* OUTPUT ARGUMENTS */
1808 /* --------------------- */
1809 /* VECOUT : VECTOR PRODUCT */
1811 /* IERCOD : ERROR CODE */
1814 /* COMMONS USED : */
1815 /* ------------------ */
1818 /* REFERENCES CALLED : */
1819 /* --------------------- */
1822 /* DESCRIPTION/NOTES/LIMITATIONS : */
1823 /* ----------------------------------- */
1825 /* ***********************************************************************
1828 /* ***********************************************************************
1833 /* ***********************************************************************
1835 /* INITIALISATIONS */
1836 /* ***********************************************************************
1839 /* Parameter adjustments */
1846 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1848 AdvApp2Var_SysBase::mgenmsg_("MMATVEC", 7L);
1852 /* ***********************************************************************
1855 /* ***********************************************************************
1857 AdvApp2Var_SysBase::mvriraz_((integer *)nligne,
1858 (char *)&vecout[1]);
1860 for (i__ = *deblig; i__ <= i__1; ++i__) {
1862 jmin = gposit[i__ * 3 + 3];
1863 jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
1864 aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
1866 for (j = jmin; j <= i__2; ++j) {
1868 somme += gmatri[k] * vecin[j];
1870 vecout[i__] = somme;
1879 /* ***********************************************************************
1881 /* ERROR PROCESSING */
1882 /* ***********************************************************************
1888 /* ***********************************************************************
1890 /* RETURN CALLING PROGRAM */
1891 /* ***********************************************************************
1896 /* ___ DESALLOCATION, ... */
1898 AdvApp2Var_SysBase::maermsg_("MMATVEC", iercod, 7L);
1900 AdvApp2Var_SysBase::mgsomsg_("MMATVEC", 7L);
1906 //=======================================================================
1907 //function : mmbulld_
1909 //=======================================================================
1910 int AdvApp2Var_MathBase::mmbulld_(integer *nbcoln,
1916 /* System generated locals */
1917 integer dtabtr_dim1, dtabtr_offset, i__1, i__2;
1919 /* Local variables */
1920 static logical ldbg;
1921 static doublereal daux;
1922 static integer nite1, nite2, nchan, i1, i2;
1924 /* ***********************************************************************
1929 /* Parsing of columns of a table of integers in increasing order */
1932 /* POINT-ENTRY, PARSING */
1933 /* INPUT ARGUMENTS : */
1934 /* -------------------- */
1935 /* - NBCOLN : Number of columns in the table */
1936 /* - NBLIGN : Number of lines in the table */
1937 /* - DTABTR : Table of integers to be parsed */
1938 /* - NUMCLE : Position of the key on the column */
1940 /* OUTPUT ARGUMENTS : */
1941 /* --------------------- */
1942 /* - DTABTR : Parsed table */
1944 /* COMMONS USED : */
1945 /* ------------------ */
1948 /* REFERENCES CALLED : */
1949 /* --------------------- */
1952 /* DESCRIPTION/NOTES/LIMITATIONS : */
1953 /* ----------------------------------- */
1954 /* Particularly performant if the table is almost parsed */
1955 /* In the opposite case it is better to use MVSHELD */
1956 /* ***********************************************************************
1959 /* Parameter adjustments */
1960 dtabtr_dim1 = *nblign;
1961 dtabtr_offset = dtabtr_dim1 + 1;
1962 dtabtr -= dtabtr_offset;
1965 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1967 AdvApp2Var_SysBase::mgenmsg_("MMBULLD", 7L);
1973 /* ***********************************************************************
1976 /* ***********************************************************************
1979 /* ---->ALGORITHM in N^2 / 2 additional iteration */
1983 /* ----> Parsing from left to the right */
1987 for (i1 = nite2; i1 <= i__1; ++i1) {
1988 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1 - 1)
1991 for (i2 = 1; i2 <= i__2; ++i2) {
1992 daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
1993 dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
1995 dtabtr[i2 + i1 * dtabtr_dim1] = daux;
2004 /* ----> Parsing from right to the left */
2009 for (i1 = nite1; i1 >= i__1; --i1) {
2010 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1
2011 - 1) * dtabtr_dim1]) {
2013 for (i2 = 1; i2 <= i__2; ++i2) {
2014 daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
2015 dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
2017 dtabtr[i2 + i1 * dtabtr_dim1] = daux;
2031 /* ***********************************************************************
2033 /* ERROR PROCESSING */
2034 /* ***********************************************************************
2037 /* ----> No errors at calling functions, only tests and loops. */
2039 /* ***********************************************************************
2041 /* RETURN CALLING PROGRAM */
2042 /* ***********************************************************************
2048 AdvApp2Var_SysBase::mgsomsg_("MMBULLD", 7L);
2055 //=======================================================================
2056 //function : AdvApp2Var_MathBase::mmcdriv_
2058 //=======================================================================
2059 int AdvApp2Var_MathBase::mmcdriv_(integer *ndimen,
2068 /* System generated locals */
2069 integer courbe_dim1, courbe_offset, crvdrv_dim1, crvdrv_offset, i__1,
2072 /* Local variables */
2073 static integer i__, j, k;
2074 static doublereal mfactk, bid;
2077 /* ***********************************************************************
2082 /* Calculate matrix of a derivate curve of order IDERIV. */
2083 /* with input parameters other than output parameters. */
2088 /* COEFFICIENTS,CURVE,DERIVATE I-EME. */
2090 /* INPUT ARGUMENTS : */
2091 /* ------------------ */
2092 /* NDIMEN : Space dimension (2 or 3 in general) */
2093 /* NCOEFF : Degree +1 of the curve. */
2094 /* COURBE : Table of coefficients of the curve. */
2095 /* IDERIV : Required order of derivation : 1=1st derivate, etc... */
2097 /* OUTPUT ARGUMENTS : */
2098 /* ------------------- */
2099 /* NCOFDV : Degree +1 of the derivative of order IDERIV of the curve. */
2100 /* CRVDRV : Table of coefficients of the derivative of order IDERIV */
2103 /* COMMONS USED : */
2104 /* ---------------- */
2106 /* REFERENCES CALLED : */
2107 /* ----------------------- */
2109 /* DESCRIPTION/NOTES/LIMITATIONS : */
2110 /* ----------------------------------- */
2112 /* ---> It is possible to take as output argument the curve */
2113 /* and the number of coeff passed at input by making : */
2114 /* CALL MMCDRIV(NDIMEN,NCOEFF,COURBE,IDERIV,NCOEFF,COURBE). */
2115 /* After this call, NCOEFF does the number of coeff of the derived */
2116 /* curve the coefficients which of are stored in CURVE. */
2117 /* Attention to the coefficients of CURVE of rank superior to */
2118 /* NCOEFF : they are not set to zero. */
2120 /* ---> Algorithm : */
2121 /* The code below was written basing on the following algorithm:
2124 /* Let P(t) = a1 + a2*t + ... an*t**n. Derivate of order k of P */
2125 /* (containing n-k coefficients) is calculated as follows : */
2127 /* Pk(t) = a(k+1)*CNP(k,k)*k! */
2128 /* + a(k+2)*CNP(k+1,k)*k! * t */
2132 /* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
2133 /* ***********************************************************************
2137 /* -------------- Case when the order of derivative is -------------------
2139 /* ---------------- greater than the degree of the curve ---------------------
2142 /* **********************************************************************
2147 /* Serves to provide the coefficients of binome (Pascal's triangle). */
2151 /* Binomial coeff from 0 to 60. read only . init par block data */
2153 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
2154 /* ----------------------------------- */
2155 /* Binomial coefficients form a triangular matrix. */
2156 /* This matrix is completed in table CNP by its transposition. */
2157 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
2159 /* Initialization is done by block-data MMLLL09.RES, */
2160 /* created by program MQINICNP.FOR). */
2161 /* **********************************************************************
2166 /* ***********************************************************************
2169 /* Parameter adjustments */
2170 crvdrv_dim1 = *ndimen;
2171 crvdrv_offset = crvdrv_dim1 + 1;
2172 crvdrv -= crvdrv_offset;
2173 courbe_dim1 = *ndimen;
2174 courbe_offset = courbe_dim1 + 1;
2175 courbe -= courbe_offset;
2178 if (*ideriv >= *ncoeff) {
2180 for (i__ = 1; i__ <= i__1; ++i__) {
2181 crvdrv[i__ + crvdrv_dim1] = 0.;
2187 /* **********************************************************************
2189 /* General processing */
2190 /* **********************************************************************
2192 /* --------------------- Calculation of Factorial(IDERIV) ------------------
2198 for (i__ = 2; i__ <= i__1; ++i__) {
2203 /* ------------ Calculation of coeff of the derived of order IDERIV ----------
2205 /* ---> Attention : coefficient binomial C(n,m) is represented in */
2206 /* MCCNP by CNP(N+1,M+1). */
2209 for (j = k + 1; j <= i__1; ++j) {
2210 bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
2212 for (i__ = 1; i__ <= i__2; ++i__) {
2213 crvdrv[i__ + (j - k) * crvdrv_dim1] = bid * courbe[i__ + j *
2220 *ncofdv = *ncoeff - *ideriv;
2222 /* -------------------------------- The end -----------------------------
2229 //=======================================================================
2230 //function : AdvApp2Var_MathBase::mmcglc1_
2232 //=======================================================================
2233 int AdvApp2Var_MathBase::mmcglc1_(integer *ndimax,
2246 /* System generated locals */
2247 integer courbe_dim1, courbe_offset, i__1;
2250 /* Local variables */
2251 static integer ndec;
2252 static doublereal tdeb, tfin;
2253 static integer iter;
2254 static doublereal oldso;
2255 static integer itmax;
2256 static doublereal sottc;
2257 static integer kk, ibb;
2258 static doublereal dif, pas;
2259 static doublereal som;
2262 /* ***********************************************************************
2267 /* Allows calculating the length of an arc of curve POLYNOMIAL */
2268 /* on an interval [A,B]. */
2272 /* LENGTH,CURVE,GAUSS,PRIVATE. */
2274 /* INPUT ARGUMENTS : */
2275 /* ------------------ */
2276 /* NDIMAX : Max. number of lines of tables */
2277 /* (i.e. max. nb of polynoms). */
2278 /* NDIMEN : Dimension of the space (nb of polynoms). */
2279 /* NCOEFF : Nb of coefficients of the polynom. This is degree + 1.
2281 /* COURBE(NDIMAX,NCOEFF) : Coefficients of the curve. */
2282 /* TDEBUT : Lower limit of the interval of integration for */
2283 /* length calculation. */
2284 /* TFINAL : Upper limit of the interval of integration for */
2285 /* length calculation. */
2286 /* EPSILN : REQIRED precision for length calculation. */
2288 /* OUTPUT ARGUMENTS : */
2289 /* ------------------- */
2290 /* XLONGC : Length of the arc of curve */
2291 /* ERREUR : Precision OBTAINED for the length calculation. */
2292 /* IERCOD : Error code, 0 OK, >0 Serious error. */
2293 /* = 1 Too much iterations, the best calculated resultat */
2294 /* (is almost ERROR) */
2295 /* = 2 Pb MMLONCV (no result) */
2296 /* = 3 NDIM or NCOEFF invalid (no result) */
2298 /* COMMONS USED : */
2299 /* ---------------- */
2301 /* REFERENCES CALLED : */
2302 /* ----------------------- */
2304 /* DESCRIPTION/NOTES/LIMITATIONS : */
2305 /* ----------------------------------- */
2306 /* The polynom is actually a set of polynoms with */
2307 /* coefficients arranged in a table of 2 indices, */
2308 /* each line relative to the polynom. */
2309 /* The polynom is defined by these coefficients ordered */
2310 /* by increasing power of the variable. */
2311 /* All polynoms have the same number of coefficients (the */
2314 /* This program cancels and replaces LENGCV, MLONGC and MLENCV. */
2316 /* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
2319 /* ***********************************************************************
2322 /* Name of the routine */
2325 /* ------------------------ General Initialization ---------------------
2328 /* Parameter adjustments */
2329 courbe_dim1 = *ndimax;
2330 courbe_offset = courbe_dim1 + 1;
2331 courbe -= courbe_offset;
2334 ibb = AdvApp2Var_SysBase::mnfndeb_();
2336 AdvApp2Var_SysBase::mgenmsg_("MMCGLC1", 7L);
2343 /* ------ Test of equity of limits */
2345 if (*tdebut == *tfinal) {
2350 /* ------ Test of the dimension and the number of coefficients */
2352 if (*ndimen <= 0 || *ncoeff <= 0) {
2356 /* ----- Nb of current cutting, nb of iteration, */
2357 /* max nb of iterations */
2364 /* ------ Variation of the nb of intervals */
2365 /* Multiplied by 2 at each iteration */
2368 pas = (*tfinal - *tdebut) / ndec;
2371 /* ------ Loop on all current NDEC intervals */
2374 for (kk = 1; kk <= i__1; ++kk) {
2376 /* ------ Limits of the current integration interval */
2378 tdeb = *tdebut + (kk - 1) * pas;
2380 mmloncv_(ndimax, ndimen, ncoeff, &courbe[courbe_offset], &tdeb, &tfin,
2392 /* ----------------- Test of the maximum number of iterations ------------
2395 /* Test if passes at least once ** */
2404 /* ------ Take into account DIF - Test of convergence */
2407 dif = (d__1 = sottc - oldso, advapp_abs(d__1));
2409 /* ------ If DIF is OK, leave..., otherwise: */
2411 if (dif > *epsiln) {
2413 /* ------ If nb iteration exceeded, leave */
2420 /* ------ Otherwise continue by cutting the initial interval.
2430 /* ------------------------------ THE END -------------------------------
2438 /* ---> PB in MMLONCV */
2444 /* ---> NCOEFF or NDIM invalid. */
2452 AdvApp2Var_SysBase::maermsg_("MMCGLC1", iercod, 7L);
2455 AdvApp2Var_SysBase::mgsomsg_("MMCGLC1", 7L);
2460 //=======================================================================
2461 //function : mmchole_
2463 //=======================================================================
2464 int mmchole_(integer *,//mxcoef,
2473 /* System generated locals */
2474 integer i__1, i__2, i__3;
2477 /* Builtin functions */
2480 /* Local variables */
2481 static logical ldbg;
2482 static integer kmin, i__, j, k;
2483 static doublereal somme;
2484 static integer ptini, ptcou;
2487 /* ***********************************************************************
2492 /* Produce decomposition of choleski of matrix A in S.S */
2493 /* Calculate inferior triangular matrix S. */
2497 /* RESOLUTION, MFACTORISATION, MATRIX_PROFILE, CHOLESKI */
2499 /* INPUT ARGUMENTS : */
2500 /* -------------------- */
2501 /* MXCOEF : Max number of terms in the hessian profile */
2502 /* DIMENS : Dimension of the problem */
2503 /* AMATRI(MXCOEF) : Coefficients of the matrix profile */
2504 /* APOSIT(1,*) : Distance diagonal-left extremity of the line
2506 /* APOSIT(2,*) : Position of diagonal terms in HESSIE */
2507 /* POSUIV(MXCOEF) : first line inferior not out of profile */
2509 /* OUTPUT ARGUMENTS : */
2510 /* --------------------- */
2511 /* CHOMAT(MXCOEF) : Inferior triangular matrix preserving the */
2512 /* profile of AMATRI. */
2513 /* IERCOD : error code */
2515 /* = 1 : non-defined positive matrix */
2517 /* COMMONS USED : */
2518 /* ------------------ */
2522 /* REFERENCES CALLED : */
2523 /* ---------------------- */
2525 /* DESCRIPTION/NOTES/LIMITATIONS : */
2526 /* ----------------------------------- */
2527 /* DEBUG LEVEL = 4 */
2528 /* ***********************************************************************
2531 /* ***********************************************************************
2536 /* ***********************************************************************
2538 /* INITIALISATIONS */
2539 /* ***********************************************************************
2542 /* Parameter adjustments */
2549 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
2551 AdvApp2Var_SysBase::mgenmsg_("MMCHOLE", 7L);
2555 /* ***********************************************************************
2558 /* ***********************************************************************
2562 for (j = 1; j <= i__1; ++j) {
2564 ptini = aposit[(j << 1) + 2];
2568 for (k = ptini - aposit[(j << 1) + 1]; k <= i__2; ++k) {
2569 /* Computing 2nd power */
2571 somme += d__1 * d__1;
2574 if (amatri[ptini] - somme < 1e-32) {
2577 chomat[ptini] = sqrt(amatri[ptini] - somme);
2581 while(posuiv[ptcou] > 0) {
2583 i__ = posuiv[ptcou];
2584 ptcou = aposit[(i__ << 1) + 2] - (i__ - j);
2586 /* Calculate the sum of S .S for k =1 a j-1 */
2590 i__2 = i__ - aposit[(i__ << 1) + 1], i__3 = j - aposit[(j << 1) +
2592 kmin = advapp_max(i__2,i__3);
2594 for (k = kmin; k <= i__2; ++k) {
2595 somme += chomat[aposit[(i__ << 1) + 2] - (i__ - k)] * chomat[
2596 aposit[(j << 1) + 2] - (j - k)];
2599 chomat[ptcou] = (amatri[ptcou] - somme) / chomat[ptini];
2605 /* ***********************************************************************
2607 /* ERROR PROCESSING */
2608 /* ***********************************************************************
2615 /* ***********************************************************************
2617 /* RETURN CALLING PROGRAM */
2618 /* ***********************************************************************
2623 AdvApp2Var_SysBase::maermsg_("MMCHOLE", iercod, 7L);
2625 AdvApp2Var_SysBase::mgsomsg_("MMCHOLE", 7L);
2631 //=======================================================================
2632 //function : AdvApp2Var_MathBase::mmcvctx_
2634 //=======================================================================
2635 int AdvApp2Var_MathBase::mmcvctx_(integer *ndimen,
2645 /* System generated locals */
2646 integer ctrtes_dim1, ctrtes_offset, crvres_dim1, crvres_offset,
2647 xmatri_dim1, xmatri_offset, tabaux_dim1, tabaux_offset, i__1,
2650 /* Local variables */
2651 static integer moup1, nordr;
2653 static integer ibb, ncf, ndv;
2654 static doublereal eps1;
2657 /* ***********************************************************************
2662 /* Calculate a polynomial curve checking the */
2663 /* passage constraints (interpolation) */
2664 /* from first derivatives, etc... to extremities. */
2665 /* Parameters at the extremities are supposed to be -1 and 1. */
2669 /* ALL, AB_SPECIFI::CONSTRAINTS&,INTERPOLATION,&CURVE */
2671 /* INPUT ARGUMENTS : */
2672 /* ------------------ */
2673 /* NDIMEN : Space Dimension. */
2674 /* NCOFMX : Nb of coeff. of curve CRVRES on each */
2676 /* NDERIV : Order of constraint with derivatives : */
2677 /* 0 --> interpolation simple. */
2678 /* 1 --> interpolation+constraints with 1st. */
2679 /* 2 --> cas (0)+ (1) + " " 2nd derivatives. */
2681 /* CTRTES : Table of constraints. */
2682 /* CTRTES(*,1,*) = contraints at -1. */
2683 /* CTRTES(*,2,*) = contraints at 1. */
2685 /* OUTPUT ARGUMENTS : */
2686 /* ------------------- */
2687 /* CRVRES : Resulting curve defined on (-1,1). */
2688 /* TABAUX : Auxilliary matrix. */
2689 /* XMATRI : Auxilliary matrix. */
2691 /* COMMONS UTILISES : */
2692 /* ---------------- */
2696 /* REFERENCES CALLED : */
2697 /* ---------------------- */
2699 /* MAERMSG R*8 DFLOAT MGENMSG */
2700 /* MGSOMSG MMEPS1 MMRSLW */
2703 /* DESCRIPTION/NOTES/LIMITATIONS : */
2704 /* ----------------------------------- */
2705 /* The polynom (or the curve) is calculated by solving a */
2706 /* system of linear equations. If the imposed degree is great */
2707 /* it is preferable to call a routine based on */
2708 /* Lagrange or Hermite interpolation depending on the case. */
2709 /* (for a high degree the matrix of the system can be badly */
2710 /* conditionned). */
2711 /* This routine returns a curve defined in (-1,1). */
2712 /* In general case, it is necessary to use MCVCTG. */
2714 /* ***********************************************************************
2717 /* Name of the routine */
2720 /* Parameter adjustments */
2721 crvres_dim1 = *ncofmx;
2722 crvres_offset = crvres_dim1 + 1;
2723 crvres -= crvres_offset;
2724 xmatri_dim1 = *nderiv + 1;
2725 xmatri_offset = xmatri_dim1 + 1;
2726 xmatri -= xmatri_offset;
2727 tabaux_dim1 = *nderiv + 1 + *ndimen;
2728 tabaux_offset = tabaux_dim1 + 1;
2729 tabaux -= tabaux_offset;
2730 ctrtes_dim1 = *ndimen;
2731 ctrtes_offset = ctrtes_dim1 * 3 + 1;
2732 ctrtes -= ctrtes_offset;
2735 ibb = AdvApp2Var_SysBase::mnfndeb_();
2737 AdvApp2Var_SysBase::mgenmsg_("MMCVCTX", 7L);
2740 AdvApp2Var_MathBase::mmeps1_(&eps1);
2742 /* ****************** CALCULATION OF EVEN COEFFICIENTS *********************
2744 /* ------------------------- Initialization -----------------------------
2747 nordr = *nderiv + 1;
2749 for (ncf = 1; ncf <= i__1; ++ncf) {
2750 tabaux[ncf + tabaux_dim1] = 1.;
2754 /* ---------------- Calculation of terms corresponding to derivatives -------
2758 for (ndv = 2; ndv <= i__1; ++ndv) {
2760 for (ncf = 1; ncf <= i__2; ++ncf) {
2761 tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
2762 tabaux_dim1] * (doublereal) ((ncf << 1) - ndv);
2768 /* ------------------ Writing the second member -----------------------
2773 for (ndv = 1; ndv <= i__1; ++ndv) {
2775 for (nd = 1; nd <= i__2; ++nd) {
2776 tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
2777 + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2778 * ctrtes_dim1]) / 2.;
2785 /* -------------------- Resolution of the system ---------------------------
2788 mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2789 xmatri_offset], iercod);
2794 for (nd = 1; nd <= i__1; ++nd) {
2796 for (ncf = 1; ncf <= i__2; ++ncf) {
2797 crvres[(ncf << 1) - 1 + nd * crvres_dim1] = xmatri[ncf + nd *
2804 /* ***************** CALCULATION OF UNEVEN COEFFICIENTS ********************
2806 /* ------------------------- Initialization -----------------------------
2811 for (ncf = 1; ncf <= i__1; ++ncf) {
2812 tabaux[ncf + tabaux_dim1] = 1.;
2816 /* ---------------- Calculation of terms corresponding to derivatives -------
2820 for (ndv = 2; ndv <= i__1; ++ndv) {
2822 for (ncf = 1; ncf <= i__2; ++ncf) {
2823 tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) *
2824 tabaux_dim1] * (doublereal) ((ncf << 1) - ndv + 1);
2830 /* ------------------ Writing of the second member -----------------------
2835 for (ndv = 1; ndv <= i__1; ++ndv) {
2837 for (nd = 1; nd <= i__2; ++nd) {
2838 tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1)
2839 + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2840 * ctrtes_dim1]) / 2.;
2847 /* -------------------- Solution of the system ---------------------------
2850 mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2851 xmatri_offset], iercod);
2856 for (nd = 1; nd <= i__1; ++nd) {
2858 for (ncf = 1; ncf <= i__2; ++ncf) {
2859 crvres[(ncf << 1) + nd * crvres_dim1] = xmatri[ncf + nd *
2866 /* --------------------------- The end ----------------------------------
2871 AdvApp2Var_SysBase::maermsg_("MMCVCTX", iercod, 7L);
2874 AdvApp2Var_SysBase::mgsomsg_("MMCVCTX", 7L);
2880 //=======================================================================
2881 //function : AdvApp2Var_MathBase::mmcvinv_
2883 //=======================================================================
2884 int AdvApp2Var_MathBase::mmcvinv_(integer *ndimax,
2891 /* Initialized data */
2893 static char nomprg[8+1] = "MMCVINV ";
2895 /* System generated locals */
2896 integer curve_dim1, curve_offset, curveo_dim1, curveo_offset, i__1, i__2;
2898 /* Local variables */
2899 static integer i__, nd, ibb;
2902 /* ***********************************************************************
2907 /* Inversion of arguments of the final curve. */
2911 /* SMOOTHING,CURVE */
2914 /* INPUT ARGUMENTS : */
2915 /* ------------------ */
2917 /* NDIM: Space Dimension. */
2918 /* NCOEF: Degree of the polynom. */
2919 /* CURVEO: The curve before inversion. */
2921 /* OUTPUT ARGUMENTS : */
2922 /* ------------------- */
2923 /* CURVE: The curve after inversion. */
2925 /* COMMONS USED : */
2926 /* ---------------- */
2927 /* REFERENCES APPELEES : */
2928 /* ----------------------- */
2929 /* DESCRIPTION/NOTES/LIMITATIONS : */
2930 /* ----------------------------------- */
2931 /* ***********************************************************************
2934 /* The name of the routine */
2935 /* Parameter adjustments */
2936 curve_dim1 = *ndimax;
2937 curve_offset = curve_dim1 + 1;
2938 curve -= curve_offset;
2939 curveo_dim1 = *ncoef;
2940 curveo_offset = curveo_dim1 + 1;
2941 curveo -= curveo_offset;
2945 ibb = AdvApp2Var_SysBase::mnfndeb_();
2947 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
2951 for (i__ = 1; i__ <= i__1; ++i__) {
2953 for (nd = 1; nd <= i__2; ++nd) {
2954 curve[nd + i__ * curve_dim1] = curveo[i__ + nd * curveo_dim1];
2963 //=======================================================================
2964 //function : mmcvstd_
2966 //=======================================================================
2967 int mmcvstd_(integer *ncofmx,
2975 /* System generated locals */
2976 integer courbe_dim1, crvcan_dim1, crvcan_offset, i__1, i__2, i__3;
2978 /* Local variables */
2979 static integer ndeg, i__, j, j1, nd, ibb;
2980 static doublereal bid;
2983 /* ***********************************************************************
2988 /* Transform curve defined between [-1,1] into [0,1]. */
2992 /* LIMITATION,RESTRICTION,CURVE */
2994 /* INPUT ARGUMENTS : */
2995 /* ------------------ */
2996 /* NDIMAX : Dimension of the space. */
2997 /* NDIMEN : Dimension of the curve. */
2998 /* NCOEFF : Degree of the curve. */
2999 /* CRVCAN(NCOFMX,NDIMEN): The curve is defined at the interval [-1,1]. */
3001 /* OUTPUT ARGUMENTS : */
3002 /* ------------------- */
3003 /* CURVE(NDIMAX,NCOEFF): Curve defined at the interval [0,1]. */
3005 /* COMMONS USED : */
3006 /* ---------------- */
3008 /* REFERENCES CALLED : */
3009 /* ----------------------- */
3011 /* DESCRIPTION/NOTES/LIMITATIONS : */
3012 /* ----------------------------------- */
3014 /* ***********************************************************************
3017 /* Name of the program. */
3020 /* **********************************************************************
3025 /* Provides binomial coefficients (Pascal triangle). */
3029 /* Binomial coefficient from 0 to 60. read only . init by block data */
3031 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3032 /* ----------------------------------- */
3033 /* Binomial coefficients form a triangular matrix. */
3034 /* This matrix is completed in table CNP by its transposition. */
3035 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
3037 /* Initialization is done with block-data MMLLL09.RES, */
3038 /* created by the program MQINICNP.FOR. */
3040 /* **********************************************************************
3045 /* ***********************************************************************
3048 /* Parameter adjustments */
3049 courbe_dim1 = *ndimax;
3051 crvcan_dim1 = *ncofmx;
3052 crvcan_offset = crvcan_dim1;
3053 crvcan -= crvcan_offset;
3056 ibb = AdvApp2Var_SysBase::mnfndeb_();
3058 AdvApp2Var_SysBase::mgenmsg_("MMCVSTD", 7L);
3062 /* ------------------ Construction of the resulting curve ----------------
3066 for (nd = 1; nd <= i__1; ++nd) {
3068 for (j = 0; j <= i__2; ++j) {
3071 for (i__ = j; i__ <= i__3; i__ += 2) {
3072 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
3076 courbe[nd + j * courbe_dim1] = bid;
3081 for (i__ = j1; i__ <= i__3; i__ += 2) {
3082 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j
3086 courbe[nd + j * courbe_dim1] -= bid;
3092 /* ------------------- Renormalization of the CURVE -------------------------
3097 for (i__ = 0; i__ <= i__1; ++i__) {
3099 for (nd = 1; nd <= i__2; ++nd) {
3100 courbe[nd + i__ * courbe_dim1] *= bid;
3107 /* ----------------------------- The end --------------------------------
3111 AdvApp2Var_SysBase::mgsomsg_("MMCVSTD", 7L);
3116 //=======================================================================
3117 //function : AdvApp2Var_MathBase::mmdrc11_
3119 //=======================================================================
3120 int AdvApp2Var_MathBase::mmdrc11_(integer *iordre,
3125 doublereal *mfactab)
3128 /* System generated locals */
3129 integer courbe_dim1, courbe_offset, points_dim2, points_offset, i__1,
3132 /* Local variables */
3134 static integer ndeg, i__, j, ndgcb, nd, ibb;
3137 /* **********************************************************************
3142 /* Calculation of successive derivatives of equation CURVE with */
3143 /* parameters -1, 1 from order 0 to order IORDRE */
3144 /* included. The calculation is produced without knowing the coefficients of */
3145 /* derivatives of the curve. */
3149 /* POSITIONING,EXTREMITIES,CURVE,DERIVATIVE. */
3151 /* INPUT ARGUMENTS : */
3152 /* ------------------ */
3153 /* IORDRE : Maximum order of calculation of derivatives. */
3154 /* NDIMEN : Dimension of the space. */
3155 /* NCOEFF : Number of coefficients of the curve (degree+1). */
3156 /* COURBE : Table of coefficients of the curve. */
3158 /* OUTPUT ARGUMENTS : */
3159 /* ------------------- */
3160 /* POINTS : Table of values of consecutive derivatives */
3161 /* of parameters -1.D0 and 1.D0. */
3162 /* MFACTAB : Auxiliary table for calculation of factorial(I).
3165 /* COMMONS USED : */
3166 /* ---------------- */
3169 /* REFERENCES CALLED : */
3170 /* ----------------------- */
3172 /* DESCRIPTION/NOTES/LIMITATIONS : */
3173 /* ----------------------------------- */
3175 /* ---> ATTENTION, the coefficients of the curve are */
3176 /* in a reverse order. */
3178 /* ---> The algorithm of calculation of derivatives is based on */
3179 /* generalization of Horner scheme : */
3181 /* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
3184 /* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
3186 /* aj = a(j-1).x + u(k-j) */
3187 /* bj = b(j-1).x + a(j-1) */
3188 /* cj = c(j-1).x + b(j-1) */
3190 /* So : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
3192 /* The algorithm is generalized easily for calculation of */
3199 /* Reference : D. KNUTH, "The Art of Computer Programming" */
3200 /* --------- Vol. 2/Seminumerical Algorithms */
3201 /* Addison-Wesley Pub. Co. (1969) */
3202 /* pages 423-425. */
3204 /* **********************************************************************
3207 /* Name of the routine */
3209 /* Parameter adjustments */
3210 points_dim2 = *iordre + 1;
3211 points_offset = (points_dim2 << 1) + 1;
3212 points -= points_offset;
3213 courbe_dim1 = *ncoeff;
3214 courbe_offset = courbe_dim1;
3215 courbe -= courbe_offset;
3218 ibb = AdvApp2Var_SysBase::mnfndeb_();
3220 AdvApp2Var_SysBase::mgenmsg_("MMDRC11", 7L);
3223 if (*iordre < 0 || *ncoeff < 1) {
3227 /* ------------------- Initialization of table POINTS -----------------
3230 ndgcb = *ncoeff - 1;
3232 for (nd = 1; nd <= i__1; ++nd) {
3233 points[(nd * points_dim2 << 1) + 1] = courbe[ndgcb + nd * courbe_dim1]
3235 points[(nd * points_dim2 << 1) + 2] = courbe[ndgcb + nd * courbe_dim1]
3241 for (nd = 1; nd <= i__1; ++nd) {
3243 for (j = 1; j <= i__2; ++j) {
3244 points[((j + nd * points_dim2) << 1) + 1] = 0.;
3245 points[((j + nd * points_dim2) << 1) + 2] = 0.;
3251 /* Calculation with parameter -1 and 1 */
3254 for (nd = 1; nd <= i__1; ++nd) {
3256 for (ndeg = 1; ndeg <= i__2; ++ndeg) {
3257 for (i__ = *iordre; i__ >= 1; --i__) {
3258 points[((i__ + nd * points_dim2) << 1) + 1] = -points[((i__ + nd
3259 * points_dim2) << 1) + 1] + points[((i__ - 1 + nd *
3260 points_dim2) << 1) + 1];
3261 points[((i__ + nd * points_dim2) << 1) + 2] += points[((i__ - 1
3262 + nd * points_dim2) << 1) + 2];
3265 points[(nd * points_dim2 << 1) + 1] = -points[(nd * points_dim2 <<
3266 1) + 1] + courbe[ndgcb - ndeg + nd * courbe_dim1];
3267 points[(nd * points_dim2 << 1) + 2] += courbe[ndgcb - ndeg + nd *
3274 /* --------------------- Multiplication by factorial(I) --------------
3278 mfac_(&mfactab[1], iordre);
3281 for (nd = 1; nd <= i__1; ++nd) {
3283 for (i__ = 2; i__ <= i__2; ++i__) {
3284 points[((i__ + nd * points_dim2) << 1) + 1] = mfactab[i__] *
3285 points[((i__ + nd * points_dim2) << 1) + 1];
3286 points[((i__ + nd * points_dim2) << 1) + 2] = mfactab[i__] *
3287 points[((i__ + nd * points_dim2) << 1) + 2];
3294 /* ---------------------------- End -------------------------------------
3299 AdvApp2Var_SysBase::mgsomsg_("MMDRC11", 7L);
3304 //=======================================================================
3305 //function : mmdrvcb_
3307 //=======================================================================
3308 int mmdrvcb_(integer *ideriv,
3317 /* System generated locals */
3318 integer courbe_dim1, tabpnt_dim1, i__1, i__2, i__3;
3320 /* Local variables */
3321 static integer ndeg, i__, j, nd, ndgcrb, iptpnt, ibb;
3324 /* ***********************************************************************
3328 /* Calculation of successive derivatives of equation CURVE with */
3329 /* parameter TPARAM from order 0 to order IDERIV included. */
3330 /* The calculation is produced without knowing the coefficients of */
3331 /* derivatives of the CURVE. */
3335 /* POSITIONING,PARAMETER,CURVE,DERIVATIVE. */
3337 /* INPUT ARGUMENTS : */
3338 /* ------------------ */
3339 /* IORDRE : Maximum order of calculation of derivatives. */
3340 /* NDIMEN : Dimension of the space. */
3341 /* NCOEFF : Number of coefficients of the curve (degree+1). */
3342 /* COURBE : Table of coefficients of the curve. */
3343 /* TPARAM : Value of the parameter where the curve should be evaluated. */
3345 /* OUTPUT ARGUMENTS : */
3346 /* ------------------- */
3347 /* TABPNT : Table of values of consecutive derivatives */
3348 /* of parameter TPARAM. */
3349 /* IERCOD : 0 = OK, */
3350 /* 1 = incoherent input. */
3352 /* COMMONS USED : */
3353 /* ---------------- */
3356 /* REFERENCES CALLED : */
3357 /* ----------------------- */
3359 /* DESCRIPTION/NOTES/LIMITATIONS : */
3360 /* ----------------------------------- */
3362 /* The algorithm of calculation of derivatives is based on */
3363 /* generalization of the Horner scheme : */
3365 /* Let C(t) = uk.t + ... + u2.t + u1.t + u0 . */
3368 /* a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
3370 /* aj = a(j-1).x + u(k-j) */
3371 /* bj = b(j-1).x + a(j-1) */
3372 /* cj = c(j-1).x + b(j-1) */
3374 /* So, it is obtained : C(x) = ak, C'(x) = bk, C"(x) = 2.ck . */
3376 /* The algorithm can be easily generalized for the calculation of */
3383 /* Reference : D. KNUTH, "The Art of Computer Programming" */
3384 /* --------- Vol. 2/Seminumerical Algorithms */
3385 /* Addison-Wesley Pub. Co. (1969) */
3386 /* pages 423-425. */
3388 /* ---> To evaluare derivatives at 0 and 1, it is preferable */
3389 /* to use routine MDRV01.FOR . */
3391 /* **********************************************************************
3394 /* Name of the routine */
3396 /* Parameter adjustments */
3397 tabpnt_dim1 = *ndim;
3399 courbe_dim1 = *ndim;
3403 ibb = AdvApp2Var_SysBase::mnfndeb_();
3405 AdvApp2Var_SysBase::mgenmsg_("MMDRVCB", 7L);
3408 if (*ideriv < 0 || *ncoeff < 1) {
3414 /* ------------------- Initialization of table TABPNT -----------------
3417 ndgcrb = *ncoeff - 1;
3419 for (nd = 1; nd <= i__1; ++nd) {
3420 tabpnt[nd] = courbe[nd + ndgcrb * courbe_dim1];
3427 iptpnt = *ndim * *ideriv;
3428 AdvApp2Var_SysBase::mvriraz_((integer *)&iptpnt,
3429 (char *)&tabpnt[tabpnt_dim1 + 1]);
3432 /* ------------------------ Calculation of parameter TPARAM ------------------
3436 for (ndeg = 1; ndeg <= i__1; ++ndeg) {
3438 for (nd = 1; nd <= i__2; ++nd) {
3439 for (i__ = *ideriv; i__ >= 1; --i__) {
3440 tabpnt[nd + i__ * tabpnt_dim1] = tabpnt[nd + i__ *
3441 tabpnt_dim1] * *tparam + tabpnt[nd + (i__ - 1) *
3445 tabpnt[nd] = tabpnt[nd] * *tparam + courbe[nd + (ndgcrb - ndeg) *
3452 /* --------------------- Multiplication by factorial(I) -------------
3456 for (i__ = 2; i__ <= i__1; ++i__) {
3458 for (j = 2; j <= i__2; ++j) {
3460 for (nd = 1; nd <= i__3; ++nd) {
3461 tabpnt[nd + i__ * tabpnt_dim1] = (doublereal) j * tabpnt[nd +
3470 /* --------------------------- The end ---------------------------------
3475 AdvApp2Var_SysBase::maermsg_("MMDRVCB", iercod, 7L);
3480 //=======================================================================
3481 //function : AdvApp2Var_MathBase::mmdrvck_
3483 //=======================================================================
3484 int AdvApp2Var_MathBase::mmdrvck_(integer *ncoeff,
3492 /* Initialized data */
3494 static doublereal mmfack[21] = { 1.,2.,6.,24.,120.,720.,5040.,40320.,
3495 362880.,3628800.,39916800.,479001600.,6227020800.,87178291200.,
3496 1.307674368e12,2.0922789888e13,3.55687428096e14,6.402373705728e15,
3497 1.21645100408832e17,2.43290200817664e18,5.109094217170944e19 };
3499 /* System generated locals */
3500 integer courbe_dim1, courbe_offset, i__1, i__2;
3502 /* Local variables */
3503 static integer i__, j, k, nd;
3504 static doublereal mfactk, bid;
3507 /* IMPLICIT INTEGER (I-N) */
3508 /* IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
3511 /* ***********************************************************************
3516 /* Calculate the value of a derived curve of order IDERIV in */
3517 /* a point of parameter TPARAM. */
3521 /* POSITIONING,CURVE,DERIVATIVE of ORDER K. */
3523 /* INPUT ARGUMENTS : */
3524 /* ------------------ */
3525 /* NCOEFF : Degree +1 of the curve. */
3526 /* NDIMEN : Dimension of the space (2 or 3 in general) */
3527 /* COURBE : Table of coefficients of the curve. */
3528 /* IDERIV : Required order of derivation : 1=1st derivative, etc... */
3529 /* TPARAM : Value of parameter of the curve. */
3531 /* OUTPUT ARGUMENTS : */
3532 /* ------------------- */
3533 /* PNTCRB : Point of parameter TPARAM on the derivative of order */
3534 /* IDERIV of CURVE. */
3536 /* COMMONS USED : */
3537 /* ---------------- */
3540 /* REFERENCES CALLED : */
3541 /* ---------------------- */
3543 /* DESCRIPTION/NOTES/LIMITATIONS : */
3544 /* ----------------------------------- */
3546 /* The code below was written basing on the following algorithm :
3549 /* Let P(t) = a1 + a2*t + ... an*t**n. The derivative of order k of P */
3550 /* (containing n-k coefficients) is calculated as follows : */
3552 /* Pk(t) = a(k+1)*CNP(k,k)*k! */
3553 /* + a(k+2)*CNP(k+1,k)*k! * t */
3557 /* + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
3559 /* Evaluation is produced following the classic Horner scheme. */
3561 /* ***********************************************************************
3565 /* Factorials (1 to 21) caculated on VAX in R*16 */
3568 /* **********************************************************************
3573 /* Serves to provide binomial coefficients (Pascal triangle). */
3577 /* Binomial Coeff from 0 to 60. read only . init by block data */
3579 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3580 /* ----------------------------------- */
3581 /* Binomial coefficients form a triangular matrix. */
3582 /* This matrix is completed in table CNP by its transposition. */
3583 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
3585 /* Initialization is done by block-data MMLLL09.RES, */
3586 /* created by program MQINICNP.FOR. */
3588 /* **********************************************************************
3593 /* ***********************************************************************
3596 /* Parameter adjustments */
3598 courbe_dim1 = *ndimen;
3599 courbe_offset = courbe_dim1 + 1;
3600 courbe -= courbe_offset;
3604 /* -------------- Case when the order of derivative is greater than -------------------
3606 /* ---------------- the degree of the curve ---------------------
3609 if (*ideriv >= *ncoeff) {
3611 for (nd = 1; nd <= i__1; ++nd) {
3617 /* **********************************************************************
3619 /* General processing*/
3620 /* **********************************************************************
3622 /* --------------------- Calculation of Factorial(IDERIV) ------------------
3626 if (*ideriv <= 21 && *ideriv > 0) {
3627 mfactk = mmfack[k - 1];
3631 for (i__ = 2; i__ <= i__1; ++i__) {
3637 /* ------- Calculation of derivative of order IDERIV of CURVE in TPARAM -----
3639 /* ---> Attention : binomial coefficient C(n,m) is represented in */
3640 /* MCCNP by CNP(N,M). */
3643 for (nd = 1; nd <= i__1; ++nd) {
3644 pntcrb[nd] = courbe[nd + *ncoeff * courbe_dim1] * mmcmcnp_.cnp[*
3645 ncoeff - 1 + k * 61] * mfactk;
3650 for (j = *ncoeff - 1; j >= i__1; --j) {
3651 bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
3653 for (nd = 1; nd <= i__2; ++nd) {
3654 pntcrb[nd] = pntcrb[nd] * *tparam + courbe[nd + j * courbe_dim1] *
3661 /* -------------------------------- The end -----------------------------
3669 //=======================================================================
3670 //function : AdvApp2Var_MathBase::mmeps1_
3672 //=======================================================================
3673 int AdvApp2Var_MathBase::mmeps1_(doublereal *epsilo)
3676 /* ***********************************************************************
3681 /* Extraction of EPS1 from COMMON MPRCSN. EPS1 is spatial zero */
3682 /* equal to 1.D-9 */
3686 /* MPRCSN,PRECISON,EPS1. */
3688 /* INPUT ARGUMENTS : */
3689 /* ------------------ */
3692 /* OUTPUT ARGUMENTS : */
3693 /* ------------------- */
3694 /* EPSILO : Value of EPS1 (spatial zero (10**-9)) */
3696 /* COMMONS USED : */
3697 /* ---------------- */
3699 /* REFERENCES CALLED : */
3700 /* ----------------------- */
3702 /* DESCRIPTION/NOTES/LIMITATIONS : */
3703 /* ----------------------------------- */
3704 /* EPS1 is ABSOLUTE spatial zero, so it is necessary */
3705 /* to use it whenever it is necessary to test if a variable */
3706 /* is null. For example, if the norm of a vector is lower than */
3707 /* EPS1, this vector is NULL ! (when one works in */
3708 /* REAL*8) It is absolutely not advised to test arguments */
3709 /* compared to EPS1**2. Taking into account the rounding errors inevitable */
3710 /* during calculations, this causes testing compared to 0.D0. */
3712 /* ***********************************************************************
3717 /* ***********************************************************************
3722 /* Gives tolerances of invalidity in stream */
3723 /* as well as limits of iterative processes */
3725 /* general context, modifiable by the user */
3729 /* PARAMETER , TOLERANCE */
3731 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
3732 /* ----------------------------------- */
3733 /* INITIALISATION : profile , **VIA MPRFTX** at input in stream
3734 /* loading of default values of the profile in MPRFTX at input */
3735 /* in stream. They are preserved in local variables of MPRFTX */
3737 /* Reset of default values : MDFINT */
3738 /* Interactive modification by the user : MDBINT */
3740 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
3741 /* MEPSPB ... EPS3,EPS4 */
3742 /* MEPSLN ... EPS2, NITERM , NITERR */
3743 /* MEPSNR ... EPS2 , NITERM */
3744 /* MITERR ... NITERR */
3746 /* ***********************************************************************
3749 /* NITERM : max nb of iterations */
3750 /* NITERR : nb of rapid iterations */
3751 /* EPS1 : tolerance of 3D null distance */
3752 /* EPS2 : tolerance of parametric null distance */
3753 /* EPS3 : tolerance to avoid division by 0.. */
3754 /* EPS4 : angular tolerance */
3758 /* ***********************************************************************
3760 *epsilo = mmprcsn_.eps1;
3765 //=======================================================================
3766 //function : mmexthi_
3768 //=======================================================================
3769 int mmexthi_(integer *ndegre,
3773 /* System generated locals */
3776 /* Local variables */
3777 static integer iadd, ideb, ndeg2, nmod2, ii, ibb;
3780 /* **********************************************************************
3785 /* Extract of common LDGRTL the weight of formulas of */
3786 /* Gauss quadrature on all roots of Legendre polynoms of degree */
3787 /* NDEGRE defined on [-1,1]. */
3791 /* ALL, AB_SPECIFI::COMMON&, EXTRACTION, &WEIGHT, &GAUSS. */
3793 /* INPUT ARGUMENTS : */
3794 /* ------------------ */
3795 /* NDEGRE : Mathematic degree of Legendre polynom. It should have */
3796 /* 2 <= NDEGRE <= 61. */
3798 /* OUTPUT ARGUMENTS : */
3799 /* ------------------- */
3800 /* HWGAUS : The table of weights of Gauss quadrature formulas */
3801 /* relative to NDEGRE roots of a polynome de Legendre de */
3804 /* COMMONS UTILISES : */
3805 /* ---------------- */
3808 /* REFERENCES CALLED : */
3809 /* ----------------------- */
3811 /* DESCRIPTION/NOTES/LIMITATIONS : */
3812 /* ----------------------------------- */
3813 /* ATTENTION: The condition on NDEGRE ( 2 <= NDEGRE <= 61) is not */
3814 /* tested. The caller should make the test.
3816 /* Name of the routine */
3819 /* Common MLGDRTL: */
3820 /* This common includes POSITIVE roots of Legendre polynims */
3821 /* AND weights of Gauss quadrature formulas on all */
3822 /* POSITIVE roots of Legendre polynoms. */
3826 /* ***********************************************************************
3831 /* The common of Legendre roots. */
3837 /* DESCRIPTION/NOTES/LIMITATIONS : */
3838 /* ----------------------------------- */
3840 /* ***********************************************************************
3846 /* ROOTAB : Table of all roots of Legendre polynoms */
3847 /* within the interval [0,1]. They are ranked for the degrees increasing from */
3849 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
3850 /* The adressing is the same. */
3851 /* HI0TAB : Table of Legendre interpolators for root x=0 */
3852 /* of polynoms of UNEVEN degree. */
3853 /* RTLTB0 : Table of Li(uk) where uk are the roots of */
3854 /* Legendre polynom of EVEN degree. */
3855 /* RTLTB1 : Table of Li(uk) where uk are the roots of */
3856 /* Legendre polynom of UNEVEN degree. */
3859 /************************************************************************
3861 /* Parameter adjustments */
3865 ibb = AdvApp2Var_SysBase::mnfndeb_();
3867 AdvApp2Var_SysBase::mgenmsg_("MMEXTHI", 7L);
3870 ndeg2 = *ndegre / 2;
3871 nmod2 = *ndegre % 2;
3873 /* Address of Gauss weight associated to the 1st strictly */
3874 /* positive root of Legendre polynom of degree NDEGRE in MLGDRTL. */
3876 iadd = ndeg2 * (ndeg2 - 1) / 2 + 1;
3878 /* Index of the 1st HWGAUS element associated to the 1st strictly */
3879 /* positive root of Legendre polynom of degree NDEGRE. */
3881 ideb = (*ndegre + 1) / 2 + 1;
3883 /* Reading of weights associated to strictly positive roots. */
3886 for (ii = ideb; ii <= i__1; ++ii) {
3887 kpt = iadd + ii - ideb;
3888 hwgaus[ii] = mlgdrtl_.hiltab[kpt + nmod2 * 465 - 1];
3892 /* For strictly negative roots, the weight is the same. */
3893 /* i.e HW(1) = HW(NDEGRE), HW(2) = HW(NDEGRE-1), etc... */
3896 for (ii = 1; ii <= i__1; ++ii) {
3897 hwgaus[ii] = hwgaus[*ndegre + 1 - ii];
3901 /* Case of uneven NDEGRE, 0 is root of Legendre polynom, */
3902 /* associated Gauss weights are loaded. */
3905 hwgaus[ndeg2 + 1] = mlgdrtl_.hi0tab[ndeg2];
3908 /* --------------------------- The end ----------------------------------
3912 AdvApp2Var_SysBase::mgsomsg_("MMEXTHI", 7L);
3917 //=======================================================================
3918 //function : mmextrl_
3920 //=======================================================================
3921 int mmextrl_(integer *ndegre,
3924 /* System generated locals */
3927 /* Local variables */
3928 static integer iadd, ideb, ndeg2, nmod2, ii, ibb;
3932 /* **********************************************************************
3937 /* Extract of the Common LDGRTL of Legendre polynom roots */
3938 /* of degree NDEGRE defined on [-1,1]. */
3942 /* ALL, AB_SPECIFI::COMMON&, EXTRACTION, &ROOT, &LEGENDRE. */
3944 /* INPUT ARGUMENTS : */
3945 /* ------------------ */
3946 /* NDEGRE : Mathematic degree of Legendre polynom. */
3947 /* It is required to have 2 <= NDEGRE <= 61. */
3949 /* OUTPUT ARGUMENTS : */
3950 /* ------------------- */
3951 /* ROOTLG : The table of roots of Legendre polynom of degree */
3952 /* NDEGRE defined on [-1,1]. */
3954 /* COMMONS USED : */
3955 /* ---------------- */
3958 /* REFERENCES CALLED : */
3959 /* ----------------------- */
3961 /* DESCRIPTION/NOTES/LIMITATIONS : */
3962 /* ----------------------------------- */
3963 /* ATTENTION: Condition of NDEGRE ( 2 <= NDEGRE <= 61) is not */
3964 /* tested. The caller should make the test. */
3966 /* **********************************************************************
3970 /* Name of the routine */
3973 /* Common MLGDRTL: */
3974 /* This common includes POSITIVE roots of Legendre polynoms */
3975 /* AND the weight of Gauss quadrature formulas on all */
3976 /* POSITIVE roots of Legendre polynoms. */
3978 /* ***********************************************************************
3983 /* The common of Legendre roots. */
3990 /* ***********************************************************************
3993 /* ROOTAB : Table of all roots of Legendre polynoms */
3994 /* within the interval [0,1]. They are ranked for the degrees increasing from */
3996 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
3997 /* The adressing is the same. */
3998 /* HI0TAB : Table of Legendre interpolators for root x=0 */
3999 /* of polynoms of UNEVEN degree. */
4000 /* RTLTB0 : Table of Li(uk) where uk are the roots of */
4001 /* Legendre polynom of EVEN degree. */
4002 /* RTLTB1 : Table of Li(uk) where uk are the roots of */
4003 /* Legendre polynom of UNEVEN degree. */
4006 /************************************************************************
4008 /* Parameter adjustments */
4012 ibb = AdvApp2Var_SysBase::mnfndeb_();
4014 AdvApp2Var_SysBase::mgenmsg_("MMEXTRL", 7L);
4017 ndeg2 = *ndegre / 2;
4018 nmod2 = *ndegre % 2;
4020 /* Address of the 1st strictly positive root of Legendre polynom */
4021 /* of degree NDEGRE in MLGDRTL. */
4023 iadd = ndeg2 * (ndeg2 - 1) / 2 + 1;
4025 /* Indice, in ROOTLG, of the 1st strictly positive root */
4026 /* of Legendre polynom of degree NDEGRE. */
4028 ideb = (*ndegre + 1) / 2 + 1;
4030 /* Reading of strictly positive roots. */
4033 for (ii = ideb; ii <= i__1; ++ii) {
4034 kpt = iadd + ii - ideb;
4035 rootlg[ii] = mlgdrtl_.rootab[kpt + nmod2 * 465 - 1];
4039 /* Strictly negative roots are equal to positive roots
4041 /* to the sign i.e RT(1) = -RT(NDEGRE), RT(2) = -RT(NDEGRE-1), etc...
4045 for (ii = 1; ii <= i__1; ++ii) {
4046 rootlg[ii] = -rootlg[*ndegre + 1 - ii];
4050 /* Case NDEGRE uneven, 0 is root of Legendre polynom. */
4053 rootlg[ndeg2 + 1] = 0.;
4056 /* -------------------------------- THE END -----------------------------
4060 AdvApp2Var_SysBase::mgenmsg_("MMEXTRL", 7L);
4065 //=======================================================================
4066 //function : AdvApp2Var_MathBase::mmfmca8_
4068 //=======================================================================
4069 int AdvApp2Var_MathBase::mmfmca8_(integer *ndimen,
4079 /* System generated locals */
4080 integer tabini_dim1, tabini_dim2, tabini_offset, tabres_dim1, tabres_dim2,
4083 /* Local variables */
4084 static integer i__, j, k, ilong;
4088 /* **********************************************************************
4093 /* Expansion of a table containing only most important things into a */
4094 /* greater data table. */
4098 /* ALL, MATH_ACCES:: CARREAU&, DECOMPRESSION, &CARREAU */
4100 /* INPUT ARGUMENTS : */
4101 /* ------------------ */
4102 /* NDIMEN: Dimension of the workspace. */
4103 /* NCOEFU: Degree +1 of the table by u. */
4104 /* NCOEFV: Degree +1 of the table by v. */
4105 /* NDIMAX: Max dimension of the space. */
4106 /* NCFUMX: Max Degree +1 of the table by u. */
4107 /* NCFVMX: Max Degree +1 of the table by v. */
4108 /* TABINI: The table to be decompressed. */
4110 /* OUTPUT ARGUMENTS : */
4111 /* ------------------- */
4112 /* TABRES: Decompressed table. */
4114 /* COMMONS USED : */
4115 /* ---------------- */
4117 /* REFERENCES CALLED : */
4118 /* ----------------------- */
4120 /* DESCRIPTION/NOTES/LIMITATIONS : */
4121 /* ----------------------------------- */
4122 /* The following call : */
4124 /* CALL MMFMCA8(NDIMEN,NCOEFU,NCOEFV,NDIMAX,NCFUMX,NCFVMX,TABINI,TABINI)
4127 /* where TABINI is input/output argument, is possible provided */
4128 /* that the caller has declared TABINI in (NDIMAX,NCFUMX,NCFVMX) */
4130 /* ATTENTION : it is not checked that NDIMAX >= NDIMEN, */
4131 /* NCOEFU >= NCFMXU and NCOEFV >= NCFMXV. */
4133 /* **********************************************************************
4137 /* Parameter adjustments */
4138 tabini_dim1 = *ndimen;
4139 tabini_dim2 = *ncoefu;
4140 tabini_offset = tabini_dim1 * (tabini_dim2 + 1) + 1;
4141 tabini -= tabini_offset;
4142 tabres_dim1 = *ndimax;
4143 tabres_dim2 = *ncfumx;
4144 tabres_offset = tabres_dim1 * (tabres_dim2 + 1) + 1;
4145 tabres -= tabres_offset;
4148 if (*ndimax == *ndimen) {
4152 /* ----------------------- decompression NDIMAX<>NDIMEN -----------------
4155 for (k = *ncoefv; k >= 1; --k) {
4156 for (j = *ncoefu; j >= 1; --j) {
4157 for (i__ = *ndimen; i__ >= 1; --i__) {
4158 tabres[i__ + (j + k * tabres_dim2) * tabres_dim1] = tabini[
4159 i__ + (j + k * tabini_dim2) * tabini_dim1];
4168 /* ----------------------- decompression NDIMAX=NDIMEN ------------------
4172 if (*ncoefu == *ncfumx) {
4175 ilong = (*ndimen << 3) * *ncoefu;
4176 for (k = *ncoefv; k >= 1; --k) {
4177 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
4178 (char *)&tabini[(k * tabini_dim2 + 1) * tabini_dim1 + 1],
4179 (char *)&tabres[(k * tabres_dim2 + 1) * tabres_dim1 + 1]);
4184 /* ----------------- decompression NDIMAX=NDIMEN,NCOEFU=NCFUMX ----------
4188 ilong = (*ndimen << 3) * *ncoefu * *ncoefv;
4189 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
4190 (char *)&tabini[tabini_offset],
4191 (char *)&tabres[tabres_offset]);
4194 /* ---------------------------- The end ---------------------------------
4201 //=======================================================================
4202 //function : AdvApp2Var_MathBase::mmfmca9_
4204 //=======================================================================
4205 int AdvApp2Var_MathBase::mmfmca9_(integer *ndimax,
4215 /* System generated locals */
4216 integer tabini_dim1, tabini_dim2, tabini_offset, tabres_dim1, tabres_dim2,
4217 tabres_offset, i__1, i__2, i__3;
4219 /* Local variables */
4220 static integer i__, j, k, ilong;
4224 /* **********************************************************************
4229 /* Compression of a data table in a table */
4230 /* containing only the main data (the input table is not removed). */
4234 /* ALL, MATH_ACCES:: CARREAU&, COMPRESSION, &CARREAU */
4236 /* INPUT ARGUMENTS : */
4237 /* ------------------ */
4238 /* NDIMAX: Max dimension of the space. */
4239 /* NCFUMX: Max degree +1 of the table by u. */
4240 /* NCFVMX: Max degree +1 of the table by v. */
4241 /* NDIMEN: Dimension of the workspace. */
4242 /* NCOEFU: Degree +1 of the table by u. */
4243 /* NCOEFV: Degree +1 of the table by v. */
4244 /* TABINI: The table to compress. */
4246 /* OUTPUT ARGUMENTS : */
4247 /* ------------------- */
4248 /* TABRES: The compressed table. */
4250 /* COMMONS USED : */
4251 /* ---------------- */
4253 /* REFERENCES CALLED : */
4254 /* ----------------------- */
4256 /* DESCRIPTION/NOTES/LIMITATIONS : */
4257 /* ----------------------------------- */
4258 /* The following call : */
4260 /* CALL MMFMCA9(NDIMAX,NCFUMX,NCFVMX,NDIMEN,NCOEFU,NCOEFV,TABINI,TABINI)
4263 /* where TABINI is input/output argument, is possible provided */
4264 /* that the caller has checked that : */
4266 /* NDIMAX > NDIMEN, */
4267 /* or NDIMAX = NDIMEN and NCFUMX > NCOEFU */
4268 /* or NDIMAX = NDIMEN, NCFUMX = NCOEFU and NCFVMX > NCOEFV */
4270 /* These conditions are not tested in the program. */
4273 /* **********************************************************************
4277 /* Parameter adjustments */
4278 tabini_dim1 = *ndimax;
4279 tabini_dim2 = *ncfumx;
4280 tabini_offset = tabini_dim1 * (tabini_dim2 + 1) + 1;
4281 tabini -= tabini_offset;
4282 tabres_dim1 = *ndimen;
4283 tabres_dim2 = *ncoefu;
4284 tabres_offset = tabres_dim1 * (tabres_dim2 + 1) + 1;
4285 tabres -= tabres_offset;
4288 if (*ndimen == *ndimax) {
4292 /* ----------------------- Compression NDIMEN<>NDIMAX -------------------
4296 for (k = 1; k <= i__1; ++k) {
4298 for (j = 1; j <= i__2; ++j) {
4300 for (i__ = 1; i__ <= i__3; ++i__) {
4301 tabres[i__ + (j + k * tabres_dim2) * tabres_dim1] = tabini[
4302 i__ + (j + k * tabini_dim2) * tabini_dim1];
4311 /* ----------------------- Compression NDIMEN=NDIMAX --------------------
4315 if (*ncoefu == *ncfumx) {
4318 ilong = (*ndimen << 3) * *ncoefu;
4320 for (k = 1; k <= i__1; ++k) {
4321 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
4322 (char *)&tabini[(k * tabini_dim2 + 1) * tabini_dim1 + 1],
4323 (char *)&tabres[(k * tabres_dim2 + 1) * tabres_dim1 + 1]);
4328 /* ----------------- Compression NDIMEN=NDIMAX,NCOEFU=NCFUMX ------------
4332 ilong = (*ndimen << 3) * *ncoefu * *ncoefv;
4333 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
4334 (char *)&tabini[tabini_offset],
4335 (char *)&tabres[tabres_offset]);
4338 /* ---------------------------- The end ---------------------------------
4345 //=======================================================================
4346 //function : AdvApp2Var_MathBase::mmfmcar_
4348 //=======================================================================
4349 int AdvApp2Var_MathBase::mmfmcar_(integer *ndimen,
4362 static integer c__8 = 8;
4363 /* System generated locals */
4364 integer patold_dim1, patold_dim2, patnew_dim1, patnew_dim2,
4365 i__1, patold_offset,patnew_offset;
4367 /* Local variables */
4368 static doublereal tbaux[1];
4369 static integer ksize, numax, kk;
4370 static long int iofst;
4371 static integer ibb, ier;
4373 /* ***********************************************************************
4378 /* LIMITATION OF A SQUARE DEFINED ON (0,1)*(0,1) BETWEEN ISOS */
4379 /* UPARA1 AND UPARA2 (BY U) AND VPARA1 AND VPARA2 BY V. */
4383 /* LIMITATION , SQUARE , PARAMETER */
4385 /* INPUT ARGUMENTS : */
4386 /* ------------------ */
4387 /* NCOFMX: MAX NUMBER OF COEFF OF THE SQUARE BY U */
4388 /* NCOEFU: NUMBER OF COEFF OF THE SQUARE BY U */
4389 /* NCOEFV: NUMBER OF COEFF OF THE SQUARE BY V */
4390 /* PATOLD : THE SQUARE IS LIMITED BY UPARA1,UPARA2 AND VPARA1,VPARA2
4392 /* UPARA1 : LOWER LIMIT OF U */
4393 /* UPARA2 : UPPER LIMIT OF U */
4394 /* VPARA1 : LOWER LIMIT OF V */
4395 /* VPARA2 : UPPER LIMIT OF V */
4397 /* OUTPUT ARGUMENTS : */
4398 /* ------------------- */
4399 /* PATNEW : RELIMITED SQUARE, DEFINED ON (0,1)**2 */
4400 /* IERCOD : =10 COEFF NB TOO GREAT OR NULL */
4401 /* =13 PB IN THE DYNAMIC ALLOCATION */
4404 /* COMMONS USED : */
4405 /* ---------------- */
4407 /* DESCRIPTION/NOTES/LIMITATIONS : */
4408 /* ----------------------------------- */
4409 /* ---> The following call : */
4410 /* CALL MMFMCAR(NCOFMX,NCOEFU,NCOEFV,PATOLD,UPARA1,UPARA2,VPARA1,VPARA2
4413 /* where PATOLD is input/output argument is absolutely legal. */
4415 /* ---> The max number of coeff by u and v of PATOLD is 61 */
4417 /* ---> If NCOEFU < NCOFMX, the data is compressed by MMFMCA9 before
4418 /* limitation by v to get time during the execution */
4419 /* of MMARC41 that follows (the square is processed as a curve of
4421 /* dimension NDIMEN*NCOEFU possessing NCOEFV coefficients). */
4423 /* ***********************************************************************
4426 /* Name of the routine */
4429 /* Parameter adjustments */
4430 patnew_dim1 = *ndimen;
4431 patnew_dim2 = *ncofmx;
4432 patnew_offset = patnew_dim1 * (patnew_dim2 + 1) + 1;
4433 patnew -= patnew_offset;
4434 patold_dim1 = *ndimen;
4435 patold_dim2 = *ncofmx;
4436 patold_offset = patold_dim1 * (patold_dim2 + 1) + 1;
4437 patold -= patold_offset;
4440 ibb = AdvApp2Var_SysBase::mnfndeb_();
4442 AdvApp2Var_SysBase::mgenmsg_("MMFMCAR", 7L);
4447 /* **********************************************************************
4449 /* TEST OF COEFFICIENT NUMBERS */
4450 /* **********************************************************************
4453 if (*ncofmx < *ncoefu) {
4457 if (*ncoefu < 1 || *ncoefu > 61 || *ncoefv < 1 || *ncoefv > 61) {
4462 /* **********************************************************************
4464 /* CASE WHEN UPARA1=VPARA1=0 AND UPARA2=VPARA2=1 */
4465 /* **********************************************************************
4468 if (*upara1 == 0. && *upara2 == 1. && *vpara1 == 0. && *vpara2 == 1.) {
4469 ksize = (*ndimen << 3) * *ncofmx * *ncoefv;
4470 AdvApp2Var_SysBase::mcrfill_((integer *)&ksize,
4471 (char *)&patold[patold_offset],
4472 (char *)&patnew[patnew_offset]);
4476 /* **********************************************************************
4478 /* LIMITATION BY U */
4479 /* **********************************************************************
4482 if (*upara1 == 0. && *upara2 == 1.) {
4486 for (kk = 1; kk <= i__1; ++kk) {
4487 mmarc41_(ndimen, ndimen, ncoefu, &patold[(kk * patold_dim2 + 1) *
4488 patold_dim1 + 1], upara1, upara2, &patnew[(kk * patnew_dim2 +
4489 1) * patnew_dim1 + 1], iercod);
4493 /* **********************************************************************
4495 /* LIMITATION BY V */
4496 /* **********************************************************************
4500 if (*vpara1 == 0. && *vpara2 == 1.) {
4504 /* ----------- LIMITATION BY V (WITH COMPRESSION I.E. NCOEFU<NCOFMX) ----
4507 numax = *ndimen * *ncoefu;
4508 if (*ncofmx != *ncoefu) {
4509 /* ------------------------- Dynamic allocation -------------------
4511 ksize = *ndimen * *ncoefu * *ncoefv;
4512 AdvApp2Var_SysBase::mcrrqst_(&c__8, &ksize, tbaux, &iofst, &ier);
4517 /* --------------- Compression by (NDIMEN,NCOEFU,NCOEFV) ------------
4519 if (*upara1 == 0. && *upara2 == 1.) {
4520 AdvApp2Var_MathBase::mmfmca9_(ndimen,
4526 &patold[patold_offset],
4529 AdvApp2Var_MathBase::mmfmca9_(ndimen,
4535 &patnew[patnew_offset],
4538 /* ------------------------- Limitation by v ------------------------
4540 mmarc41_(&numax, &numax, ncoefv, &tbaux[iofst], vpara1, vpara2, &
4541 tbaux[iofst], iercod);
4542 /* --------------------- Expansion of TBAUX into PATNEW -------------
4544 AdvApp2Var_MathBase::mmfmca8_(ndimen, ncoefu, ncoefv, ndimen, ncofmx, ncoefv, &tbaux[iofst]
4545 , &patnew[patnew_offset]);
4548 /* -------- LIMITATION BY V (WITHOUT COMPRESSION I.E. NCOEFU=NCOFMX) ---
4552 if (*upara1 == 0. && *upara2 == 1.) {
4553 mmarc41_(&numax, &numax, ncoefv, &patold[patold_offset], vpara1,
4554 vpara2, &patnew[patnew_offset], iercod);
4556 mmarc41_(&numax, &numax, ncoefv, &patnew[patnew_offset], vpara1,
4557 vpara2, &patnew[patnew_offset], iercod);
4562 /* **********************************************************************
4565 /* **********************************************************************
4570 AdvApp2Var_SysBase::mcrdelt_(&c__8, &ksize, tbaux, &iofst, &ier);
4576 /* ------------------------------ The end -------------------------------
4581 AdvApp2Var_SysBase::maermsg_("MMFMCAR", iercod, 7L);
4584 AdvApp2Var_SysBase::mgsomsg_("MMFMCAR", 7L);
4590 //=======================================================================
4591 //function : AdvApp2Var_MathBase::mmfmcb5_
4593 //=======================================================================
4594 int AdvApp2Var_MathBase::mmfmcb5_(integer *isenmsc,
4605 /* System generated locals */
4606 integer courb1_dim1, courb1_offset, courb2_dim1, courb2_offset, i__1,
4609 /* Local variables */
4610 static integer i__, nboct, nd;
4613 /* **********************************************************************
4618 /* Reformating (and eventual compression/decompression) of curve */
4619 /* (ndim,.) by (.,ndim) and vice versa. */
4623 /* ALL , MATH_ACCES :: */
4624 /* COURBE&, REORGANISATION,COMPRESSION,INVERSION , &COURBE */
4626 /* INPUT ARGUMENTS : */
4627 /* -------------------- */
4628 /* ISENMSC : required direction of the transfer : */
4629 /* 1 : passage of (NDIMEN,.) ---> (.,NDIMEN) direction to AB
4631 /* -1 : passage of (.,NDIMEN) ---> (NDIMEN,.) direction to TS,T
4633 /* NDIMAX : format / dimension */
4634 /* NCF1MX : format by t of COURB1 */
4635 /* if ISENMSC= 1 : COURB1: The curve to be processed (NDIMAX,.) */
4636 /* NCOEFF : number of coeff of the curve */
4637 /* NCF2MX : format by t of COURB2 */
4638 /* NDIMEN : dimension of the curve and format of COURB2 */
4639 /* if ISENMSC=-1 : COURB2: The curve to be processed (.,NDIMEN) */
4641 /* OUTPUT ARGUMENTS : */
4642 /* --------------------- */
4643 /* if ISENMSC= 1 : COURB2: The resulting curve (.,NDIMEN) */
4644 /* if ISENMSC=-1 : COURB1: The resulting curve (NDIMAX,.) */
4646 /* COMMONS USED : */
4647 /* ------------------ */
4649 /* REFERENCES CALLED : */
4650 /* --------------------- */
4652 /* DESCRIPTION/NOTES/LIMITATIONS : */
4653 /* ----------------------------------- */
4654 /* allow to process the usual transfers as follows : */
4655 /* | ---- ISENMSC = 1 ---- | | ---- ISENMSC =-1 ----- | */
4656 /* TS (3,21) --> (21,3) AB ; AB (21,3) --> (3,21) TS */
4657 /* TS (3,21) --> (NU,3) AB ; AB (NU,3) --> (3,21) TS */
4658 /* (3,NU) --> (21,3) AB ; AB (21,3) --> (3,NU) */
4659 /* (3,NU) --> (NU,3) AB ; AB (NU,3) --> (3,NU) */
4661 /* ***********************************************************************
4665 /* Parameter adjustments */
4666 courb1_dim1 = *ndimax;
4667 courb1_offset = courb1_dim1 + 1;
4668 courb1 -= courb1_offset;
4669 courb2_dim1 = *ncf2mx;
4670 courb2_offset = courb2_dim1 + 1;
4671 courb2 -= courb2_offset;
4674 if (*ndimen > *ndimax || *ncoeff > *ncf1mx || *ncoeff > *ncf2mx) {
4678 if (*ndimen == 1 && *ncf1mx == *ncf2mx) {
4679 nboct = *ncf2mx << 3;
4680 if (*isenmsc == 1) {
4681 AdvApp2Var_SysBase::mcrfill_((integer *)&nboct,
4682 (char *)&courb1[courb1_offset],
4683 (char *)&courb2[courb2_offset]);
4685 if (*isenmsc == -1) {
4686 AdvApp2Var_SysBase::mcrfill_((integer *)&nboct,
4687 (char *)&courb2[courb2_offset],
4688 (char *)&courb1[courb1_offset]);
4695 if (*isenmsc == 1) {
4697 for (nd = 1; nd <= i__1; ++nd) {
4699 for (i__ = 1; i__ <= i__2; ++i__) {
4700 courb2[i__ + nd * courb2_dim1] = courb1[nd + i__ *
4706 } else if (*isenmsc == -1) {
4708 for (nd = 1; nd <= i__1; ++nd) {
4710 for (i__ = 1; i__ <= i__2; ++i__) {
4711 courb1[nd + i__ * courb1_dim1] = courb2[i__ + nd *
4723 /* ***********************************************************************
4731 AdvApp2Var_SysBase::maermsg_("MMFMCB5", iercod, 7L);
4736 //=======================================================================
4737 //function : AdvApp2Var_MathBase::mmfmtb1_
4739 //=======================================================================
4740 int AdvApp2Var_MathBase::mmfmtb1_(integer *maxsz1,
4750 static integer c__8 = 8;
4752 /* System generated locals */
4753 integer table1_dim1, table1_offset, table2_dim1, table2_offset, i__1,
4756 /* Local variables */
4757 static doublereal work[1];
4758 static integer ilong, isize, ii, jj, ier;
4759 static long int iofst,iipt, jjpt;
4762 /************************************************************************
4767 /* Inversion of elements of a rectangular table (T1(i,j) */
4768 /* loaded in T2(j,i)) */
4772 /* ALL, MATH_ACCES :: TABLEAU&, INVERSION, &TABLEAU */
4774 /* INPUT ARGUMENTS : */
4775 /* ------------------ */
4776 /* MAXSZ1: Max Nb of elements by the 1st dimension of TABLE1. */
4777 /* TABLE1: Table of reals by two dimensions. */
4778 /* ISIZE1: Nb of useful elements of TABLE1 on the 1st dimension */
4779 /* JSIZE1: Nb of useful elements of TABLE1 on the 2nd dimension */
4780 /* MAXSZ2: Nb max of elements by the 1st dimension of TABLE2. */
4782 /* OUTPUT ARGUMENTS : */
4783 /* ------------------- */
4784 /* TABLE2: Table of reals by two dimensions, containing the transposition
4785 /* of the rectangular table TABLE1. */
4786 /* ISIZE2: Nb of useful elements of TABLE2 on the 1st dimension */
4787 /* JSIZE2: Nb of useful elements of TABLE2 on the 2nd dimension */
4788 /* IERCOD: Erroe coder. */
4790 /* = 1, error in the dimension of tables */
4791 /* ether MAXSZ1 < ISIZE1 (table TABLE1 too small). */
4792 /* or MAXSZ2 < JSIZE1 (table TABLE2 too small). */
4794 /* COMMONS USED : */
4795 /* ---------------- */
4797 /* REFERENCES CALLED : */
4798 /* ---------------------- */
4800 /* DESCRIPTION/NOTES/LIMITATIONS : */
4801 /* ----------------------------------- */
4802 /* It is possible to use TABLE1 as input and output table i.e. */
4804 /* CALL MMFMTB1(MAXSZ1,TABLE1,ISIZE1,JSIZE1,MAXSZ2,TABLE1 */
4805 /* ,ISIZE2,JSIZE2,IERCOD) */
4808 /* **********************************************************************
4812 /* Parameter adjustments */
4813 table1_dim1 = *maxsz1;
4814 table1_offset = table1_dim1 + 1;
4815 table1 -= table1_offset;
4816 table2_dim1 = *maxsz2;
4817 table2_offset = table2_dim1 + 1;
4818 table2 -= table2_offset;
4822 if (*isize1 > *maxsz1 || *jsize1 > *maxsz2) {
4827 isize = *maxsz2 * *isize1;
4828 AdvApp2Var_SysBase::mcrrqst_(&c__8, &isize, work, &iofst, &ier);
4833 /* DO NOT BE AFRAID OF CRUSHING. */
4836 for (ii = 1; ii <= i__1; ++ii) {
4837 iipt = (ii - 1) * *maxsz2 + iofst;
4839 for (jj = 1; jj <= i__2; ++jj) {
4840 jjpt = iipt + (jj - 1);
4841 work[jjpt] = table1[ii + jj * table1_dim1];
4847 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
4848 (char *)&work[iofst],
4849 (char *)&table2[table2_offset]);
4851 /* -------------- The number of elements of TABLE2 is returned ------------
4860 /* ------------------------------- THE END ------------------------------
4862 /* --> Invalid input. */
4866 /* --> Pb of allocation. */
4873 AdvApp2Var_SysBase::mcrdelt_(&c__8, &isize, work, &iofst, &ier);
4881 //=======================================================================
4882 //function : AdvApp2Var_MathBase::mmgaus1_
4884 //=======================================================================
4885 int AdvApp2Var_MathBase::mmgaus1_(integer *ndimf,
4902 /* System generated locals */
4905 /* Local variables */
4906 static integer ndeg;
4907 static doublereal h__[20];
4909 static doublereal t, u[20], x;
4910 static integer idimf;
4911 static doublereal c1x, c2x;
4912 /* **********************************************************************
4918 /* Calculate the integral of function BFUNX passed in parameter */
4919 /* between limits XD and XF . */
4920 /* The function should be calculated for any value */
4921 /* of the variable in the given interval.. */
4922 /* The method GAUSS-LEGENDRE is used.
4923 /* For explications refer to the book : */
4924 /* Complements de mathematiques a l'usage des Ingenieurs de */
4925 /* l'electrotechnique et des telecommunications. */
4926 /* Par Andre ANGOT - Collection technique et scientifique du CNET
4929 /* The degree of LEGENDRE polynoms used is passed in parameter.
4933 /* INTEGRATION,LEGENDRE,GAUSS */
4935 /* INPUT ARGUMENTS : */
4936 /* ------------------ */
4938 /* NDIMF : Dimension of the function */
4939 /* BFUNX : Function to integrate passed as argument */
4940 /* Should be declared as EXTERNAL in the call routine. */
4941 /* SUBROUTINE BFUNX(NDIMF,X,VAL,IER) */
4943 /* K : Parameter determining the degree of the LEGENDRE polynom that
4945 /* can take a value between 0 and 10. */
4946 /* The degree of the polynom is equal to 4 k, that is 4, 8,
4948 /* 12, 16, 20, 24, 28, 32, 36 and 40. */
4949 /* If K is not correct, the degree is set to 40 directly.
4951 /* XD : Lower limit of the interval of integration. */
4952 /* XF : Upper limit of the interval of integration. */
4953 /* SAUX1 : Auxiliary table */
4954 /* SAUX2 : Auxiliary table */
4956 /* OUTPUT ARGUMENTS : */
4957 /* ------------------- */
4959 /* SOMME : Value of the integral */
4960 /* NITER : Number of iterations to be carried out. */
4961 /* It is equal to the degree of the polynom. */
4963 /* IER : Error code : */
4964 /* < 0 ==> Attention - Warning */
4965 /* = 0 ==> Everything is OK */
4966 /* > 0 ==> Critical error - Apply special processing */
4967 /* ==> Error in the calculation of BFUNX (return code */
4968 /* of this routine */
4970 /* If error => SUM = 0 */
4972 /* COMMONS USED : */
4973 /* ----------------- */
4977 /* REFERENCES CALLED : */
4978 /* ---------------------- */
4981 /* @ BFUNX MVGAUS0 */
4983 /* DESCRIPTION/NOTES/LIMITATIONS : */
4984 /* --------------------------------- */
4986 /* See the explanations detailed in the listing */
4987 /* Use of the GAUSS method (orthogonal polynoms) */
4988 /* The symmetry of roots of these polynomes is used */
4989 /* Depending on K, the degree of the interpolated polynom grows.
4991 /* If you wish to calculate the integral with a given precision, */
4992 /* loop on k varying from 1 to 10 and test the difference of 2
4994 /* consecutive iterations. Stop the loop if this difference is less that
4995 /* an epsilon value set to 10E-6 for example. */
4996 /* If S1 and S2 are 2 successive iterations, test following this example :
4999 /* AF=DABS(S1-S2) */
5001 /* If AS < 1 test if FS < eps otherwise test if AF/AS < eps
5003 /* -- ----- ----- */
5005 /************************************************************************
5008 /************************************************************************
5013 /* ****** General Initialization */
5015 /* Parameter adjustments */
5021 AdvApp2Var_SysBase::mvriraz_((integer *)ndimf,
5025 /* ****** Loading of coefficients U and H ** */
5026 /* -------------------------------------------- */
5028 mvgaus0_(k, u, h__, &ndeg, iercod);
5033 /* ****** C1X => Medium interval point [XD,XF] */
5034 /* ****** C2X => 1/2 amplitude interval [XD,XF] */
5036 c1x = (*xf + *xd) * .5;
5037 c2x = (*xf - *xd) * .5;
5039 /* ---------------------------------------- */
5040 /* ****** Integration for degree NDEG ** */
5041 /* ---------------------------------------- */
5044 for (j = 1; j <= i__1; ++j) {
5048 (*bfunx)(ndimf, &x, &saux1[1], iercod);
5054 (*bfunx)(ndimf, &x, &saux2[1], iercod);
5060 for (idimf = 1; idimf <= i__2; ++idimf) {
5061 somme[idimf] += h__[j - 1] * (saux1[idimf] + saux2[idimf]);
5068 for (idimf = 1; idimf <= i__1; ++idimf) {
5069 somme[idimf] *= c2x;
5072 /* ****** End of sub-program ** */
5078 //=======================================================================
5079 //function : mmherm0_
5081 //=======================================================================
5082 int mmherm0_(doublereal *debfin,
5085 static integer c__576 = 576;
5086 static integer c__6 = 6;
5089 /* System generated locals */
5093 /* Local variables */
5094 static doublereal amat[36] /* was [6][6] */;
5095 static integer iord[2];
5096 static doublereal prod;
5097 static integer iord1, iord2;
5098 static doublereal miden[36] /* was [6][6] */;
5099 static integer ncmat;
5100 static doublereal epspi, d1, d2;
5101 static integer ii, jj, pp, ncf;
5102 static doublereal cof[6];
5103 static integer iof[2], ier;
5104 static doublereal mat[36] /* was [6][6] */;
5106 static doublereal abid[72] /* was [12][6] */;
5107 /* ***********************************************************************
5112 /* INIT OF COEFFS. OF POLYNOMS OF HERMIT INTERPOLATION */
5116 /* MATH_ACCES :: HERMITE */
5118 /* INPUT ARGUMENTS */
5119 /* -------------------- */
5120 /* DEBFIN : PARAMETERS DEFINING THE CONSTRAINTS */
5121 /* DEBFIN(1) : FIRST PARAMETER */
5122 /* DEBFIN(2) : SECOND PARAMETER */
5124 /* ONE SHOULD HAVE: */
5125 /* ABS (DEBFIN(I)) < 100 */
5127 /* (ABS(DEBFIN(1)+ABS(DEBFIN(2))) > 1/100 */
5128 /* (for overflows) */
5130 /* ABS(DEBFIN(2)-DEBFIN(1)) / (ABS(DEBFIN(1)+ABS(DEBFIN(2))) > 1/100
5132 /* (for the conditioning) */
5135 /* OUTPUT ARGUMENTS : */
5136 /* --------------------- */
5138 /* IERCOD : Error code : 0 : O.K. */
5139 /* 1 : value of DEBFIN */
5140 /* are unreasonable */
5141 /* -1 : init was already done */
5142 /* (OK but no processing) */
5144 /* COMMONS USED : */
5145 /* ------------------ */
5147 /* REFERENCES CALLED : */
5148 /* ---------------------- */
5151 /* DESCRIPTION/NOTES/LIMITATIONS : */
5152 /* ----------------------------------- */
5154 /* This program initializes the coefficients of Hermit polynoms */
5155 /* that are read later by MMHERM1 */
5156 /* ***********************************************************************
5161 /* **********************************************************************
5166 /* Used to STORE coefficients of Hermit interpolation polynoms
5172 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5173 /* ----------------------------------- */
5175 /* The coefficients of hermit polynoms are calculated by */
5176 /* the routine MMHERM0 and read by the routine MMHERM1 */
5178 /* **********************************************************************
5185 /* NBCOEF is the size of CMHERM (see below) */
5186 /* ***********************************************************************
5195 /* ***********************************************************************
5198 /* ***********************************************************************
5202 /* Parameter adjustments */
5206 d1 = advapp_abs(debfin[1]);
5207 if (d1 > (float)100.) {
5211 d2 = advapp_abs(debfin[2]);
5212 if (d2 > (float)100.) {
5217 if (d2 < (float).01) {
5221 d1 = (d__1 = debfin[2] - debfin[1], advapp_abs(d__1));
5222 if (d1 / d2 < (float).01) {
5227 /* ***********************************************************************
5229 /* Initialization */
5230 /* ***********************************************************************
5238 /* ***********************************************************************
5241 /* IS IT ALREADY INITIALIZED ? */
5243 d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
5246 if (debfin[1] != mmcmher_.tdebut) {
5249 if (debfin[2] != mmcmher_.tfinal) {
5252 if (d1 != mmcmher_.verifi) {
5260 /* ***********************************************************************
5263 /* ***********************************************************************
5269 /* Init. matrix identity : */
5272 AdvApp2Var_SysBase::mvriraz_((integer *)&ncmat,
5275 for (ii = 1; ii <= 6; ++ii) {
5276 miden[ii + ii * 6 - 7] = 1.;
5282 /* Init to 0 of table CMHERM */
5284 AdvApp2Var_SysBase::mvriraz_((integer *)&c__576, (char *)mmcmher_.cmherm);
5286 /* Calculation by solution of linear systems */
5288 for (iord1 = -1; iord1 <= 2; ++iord1) {
5289 for (iord2 = -1; iord2 <= 2; ++iord2) {
5296 iof[1] = iord[0] + 1;
5299 ncf = iord[0] + iord[1] + 2;
5301 /* Calculate matrix MAT to invert: */
5303 for (cot = 1; cot <= 2; ++cot) {
5306 if (iord[cot - 1] > -1) {
5309 for (jj = 1; jj <= i__1; ++jj) {
5315 i__1 = iord[cot - 1] + 1;
5316 for (pp = 1; pp <= i__1; ++pp) {
5318 ii = pp + iof[cot - 1];
5323 for (jj = 1; jj <= i__2; ++jj) {
5324 mat[ii + jj * 6 - 7] = (float)0.;
5329 for (jj = pp; jj <= i__2; ++jj) {
5331 /* everything is done in these 3 lines
5334 mat[ii + jj * 6 - 7] = cof[jj - 1] * prod;
5335 cof[jj - 1] *= jj - pp;
5336 prod *= debfin[cot];
5349 AdvApp2Var_MathBase::mmmrslwd_(&c__6, &ncf, &ncf, mat, miden, &epspi, abid, amat, &
5356 for (cot = 1; cot <= 2; ++cot) {
5357 i__1 = iord[cot - 1] + 1;
5358 for (pp = 1; pp <= i__1; ++pp) {
5360 for (ii = 1; ii <= i__2; ++ii) {
5361 mmcmher_.cmherm[ii + (pp + (cot + ((iord1 + (iord2 <<
5362 2)) << 1)) * 3) * 6 + 155] = amat[ii + (pp +
5363 iof[cot - 1]) * 6 - 7];
5376 /* ***********************************************************************
5379 /* The initialized flag is located: */
5381 mmcmher_.tdebut = debfin[1];
5382 mmcmher_.tfinal = debfin[2];
5384 d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
5385 mmcmher_.verifi = d1 * 16111959;
5388 /* ***********************************************************************
5393 /* ***********************************************************************
5404 /* ***********************************************************************
5409 AdvApp2Var_SysBase::maermsg_("MMHERM0", iercod, 7L);
5411 /* ***********************************************************************
5416 //=======================================================================
5417 //function : mmherm1_
5419 //=======================================================================
5420 int mmherm1_(doublereal *debfin,
5426 /* System generated locals */
5427 integer hermit_dim1, hermit_dim2, hermit_offset;
5429 /* Local variables */
5430 static integer nbval;
5431 static doublereal d1;
5434 /* ***********************************************************************
5439 /* reading of coeffs. of HERMIT interpolation polynoms */
5443 /* MATH_ACCES :: HERMIT */
5445 /* INPUT ARGUMENTS : */
5446 /* -------------------- */
5447 /* DEBFIN : PARAMETES DEFINING THE CONSTRAINTS */
5448 /* DEBFIN(1) : FIRST PARAMETER */
5449 /* DEBFIN(2) : SECOND PARAMETER */
5451 /* Should be equal to the corresponding arguments during the */
5452 /* last call to MMHERM0 for the initialization of coeffs. */
5454 /* ORDRMX : indicates the dimensioning of HERMIT: */
5455 /* there is no choice : ORDRMX should be equal to the value */
5456 /* of PARAMETER IORDMX of INCLUDE MMCMHER, or 2 for the moment */
5458 /* IORDRE (2) : Orders of constraints in each corresponding parameter DEBFIN(I)
5459 /* should be between -1 (no constraints) and ORDRMX. */
5462 /* OUTPUT ARGUMENTS : */
5463 /* --------------------- */
5465 /* HERMIT : HERMIT(1:IORDRE(1)+IORDRE(2)+2, j, cote) are the */
5466 /* coefficients in the canonic base of Hermit polynom */
5467 /* corresponding to orders IORDRE with parameters DEBFIN for */
5468 /* the constraint of order j on DEBFIN(cote). j is between 0 and IORDRE(cote). */
5471 /* IERCOD : Error code : */
5472 /* -1: O.K but necessary to reinitialize the coefficients */
5473 /* (info for optimization) */
5475 /* 1 : Error in MMHERM0 */
5476 /* 2 : arguments invalid */
5478 /* COMMONS USED : */
5479 /* ------------------ */
5481 /* REFERENCES CALLED : */
5482 /* ---------------------- */
5485 /* DESCRIPTION/NOTES/LIMITATIONS : */
5486 /* ----------------------------------- */
5488 /* This program reads coefficients of Hermit polynoms */
5489 /* that were earlier initialized by MMHERM0 */
5491 /* PMN : initialisation is no more done by the caller. */
5494 /* ***********************************************************************
5499 /* **********************************************************************
5504 /* Serves to STORE the coefficients of Hermit interpolation polynoms
5510 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5511 /* ----------------------------------- */
5513 /* the coefficients of Hetmit polynoms are calculated by */
5514 /* routine MMHERM0 and read by routine MMHERM1 */
5517 /* **********************************************************************
5524 /* NBCOEF is the size of CMHERM (see lower) */
5528 /* ***********************************************************************
5535 /* ***********************************************************************
5537 /* Initializations */
5538 /* ***********************************************************************
5541 /* Parameter adjustments */
5543 hermit_dim1 = (*ordrmx << 1) + 2;
5544 hermit_dim2 = *ordrmx + 1;
5545 hermit_offset = hermit_dim1 * hermit_dim2 + 1;
5546 hermit -= hermit_offset;
5553 /* ***********************************************************************
5556 /* ***********************************************************************
5564 for (cot = 1; cot <= 2; ++cot) {
5565 if (iordre[cot] < -1) {
5568 if (iordre[cot] > *ordrmx) {
5575 /* IS-IT CORRECTLY INITIALIZED ? */
5577 d1 = advapp_abs(debfin[1]) + advapp_abs(debfin[2]);
5580 /* OTHERWISE IT IS INITIALIZED */
5582 if (debfin[1] != mmcmher_.tdebut || debfin[2] != mmcmher_.tfinal || d1
5583 != mmcmher_.verifi) {
5585 mmherm0_(&debfin[1], iercod);
5592 /* ***********************************************************************
5595 /* ***********************************************************************
5600 AdvApp2Var_SysBase::msrfill_(&nbval, &mmcmher_.cmherm[((((iordre[1] + (iordre[2] << 2)) << 1)
5601 + 1) * 3 + 1) * 6 + 156], &hermit[hermit_offset]);
5603 /* ***********************************************************************
5608 /* ***********************************************************************
5619 /* ***********************************************************************
5624 AdvApp2Var_SysBase::maermsg_("MMHERM1", iercod, 7L);
5626 /* ***********************************************************************
5631 //=======================================================================
5632 //function : AdvApp2Var_MathBase::mmhjcan_
5634 //=======================================================================
5635 int AdvApp2Var_MathBase::mmhjcan_(integer *ndimen,
5646 static integer c__2 = 2;
5647 static integer c__21 = 21;
5648 /* System generated locals */
5649 integer tcbold_dim1, tcbold_dim2, tcbold_offset, tcbnew_dim1, tcbnew_dim2,
5650 tcbnew_offset, i__1, i__2, i__3, i__4, i__5;
5653 /* Local variables */
5654 static logical ldbg;
5655 static integer ndeg;
5656 static doublereal taux1[21];
5657 static integer d__, e, i__, k;
5658 static doublereal mfact;
5659 static integer ncoeff;
5660 static doublereal tjacap[21];
5661 static integer iordre[2];
5662 static doublereal hermit[36]/* was [6][3][2] */, ctenor, bornes[2];
5664 static integer aux1, aux2;
5666 /* ***********************************************************************
5671 /* CONVERSION OF TABLE TCBOLD OF POLYNOMIAL CURVE COEFFICIENTS */
5672 /* EXPRESSED IN HERMIT JACOBI BASE, INTO A */
5673 /* TABLE OF COEFFICIENTS TCBNEW OF COURVES EXPRESSED IN THE CANONIC BASE */
5677 /* CANNONIC, HERMIT, JACCOBI */
5679 /* INPUT ARGUMENTS : */
5680 /* -------------------- */
5681 /* ORDHER : ORDER OF HERMIT POLYNOMS OR ORDER OF CONTINUITY */
5682 /* NCOEFS : NUMBER OF COEFFICIENTS OF A POLYNOMIAL CURVE */
5683 /* FOR ONE OF ITS NDIM COMPONENTS;(DEGREE+1 OF THE CURVE)
5685 /* NDIM : DIMENSION OF THE CURVE */
5686 /* CBHEJA : TABLE OF COEFFICIENTS OF THE CURVE IN THE BASE */
5688 /* (H(0,-1),..,H(ORDHER,-1),H(0,1),..,H(ORDHER,1), */
5689 /* JA(ORDHER+1,2*ORDHER+2),....,JA(ORDHER+1,NCOEFS-1) */
5691 /* OUTPUT ARGUMENTS : */
5692 /* --------------------- */
5693 /* CBRCAN : TABLE OF COEFFICIENTS OF THE CURVE IN THE CANONIC BASE */
5696 /* COMMONS USED : */
5697 /* ------------------ */
5700 /* REFERENCES CALLED : */
5701 /* --------------------- */
5704 /* ***********************************************************************
5708 /* ***********************************************************************
5713 /* Providesinteger constants from 0 to 1000 */
5719 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
5720 /* ----------------------------------- */
5722 /* ***********************************************************************
5726 /* ***********************************************************************
5732 /* ***********************************************************************
5734 /* INITIALIZATION */
5735 /* ***********************************************************************
5738 /* Parameter adjustments */
5740 tcbnew_dim1 = *ndimen;
5741 tcbnew_dim2 = *ncflim;
5742 tcbnew_offset = tcbnew_dim1 * (tcbnew_dim2 + 1) + 1;
5743 tcbnew -= tcbnew_offset;
5744 tcbold_dim1 = *ndimen;
5745 tcbold_dim2 = *ncflim;
5746 tcbold_offset = tcbold_dim1 * (tcbold_dim2 + 1) + 1;
5747 tcbold -= tcbold_offset;
5750 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
5752 AdvApp2Var_SysBase::mgenmsg_("MMHJCAN", 7L);
5759 /* ***********************************************************************
5762 /* ***********************************************************************
5772 /* CALCULATION OF HERMIT POLYNOMS IN THE CANONIC BASE ON (-1,1) */
5775 iordre[0] = *orcont;
5776 iordre[1] = *orcont;
5777 mmherm1_(bornes, &c__2, iordre, hermit, &ier);
5787 for (e = 1; e <= i__1; ++e) {
5789 ctenor = (tdecop[e] - tdecop[e - 1]) / 2;
5797 for (d__ = 1; d__ <= i__2; ++d__) {
5799 /* CONVERSION OF THE COEFFICIENTS OF THE PART OF THE CURVE EXPRESSED */
5800 /* IN HERMIT BASE, INTO THE CANONIC BASE */
5802 AdvApp2Var_SysBase::mvriraz_((integer *)&ncoeff, (char *)taux1);
5805 for (k = 1; k <= i__3; ++k) {
5807 for (i__ = 1; i__ <= i__4; ++i__) {
5809 mfact = AdvApp2Var_MathBase::pow__di(&ctenor, &i__5);
5810 taux1[k - 1] += (tcbold[d__ + (i__ + e * tcbold_dim2) *
5811 tcbold_dim1] * hermit[k + (i__ + 2) * 6 - 19] +
5812 tcbold[d__ + (i__ + aux1 + e * tcbold_dim2) *
5813 tcbold_dim1] * hermit[k + (i__ + 5) * 6 - 19]) *
5820 for (i__ = aux2 + 1; i__ <= i__3; ++i__) {
5821 taux1[i__ - 1] = tcbold[d__ + (i__ + e * tcbold_dim2) *
5825 /* CONVERSION OF THE COEFFICIENTS OF THE PART OF THE CURVE EXPRESSED */
5826 /* IN CANONIC-JACOBI BASE, INTO THE CANONIC BASE */
5830 AdvApp2Var_MathBase::mmapcmp_(&minombr_.nbr[1], &c__21, &ncoeff, taux1, tjacap);
5831 AdvApp2Var_MathBase::mmjacan_(orcont, &ndeg, tjacap, taux1);
5833 /* RECOPY THE COEFS RESULTING FROM THE CONVERSION IN THE TABLE */
5837 for (i__ = 1; i__ <= i__3; ++i__) {
5838 tcbnew[d__ + (i__ + e * tcbnew_dim2) * tcbnew_dim1] = taux1[
5847 /* ***********************************************************************
5849 /* PROCESSING OF ERRORS */
5850 /* ***********************************************************************
5860 /* ***********************************************************************
5862 /* RETURN CALLING PROGRAM */
5863 /* ***********************************************************************
5868 AdvApp2Var_SysBase::maermsg_("MMHJCAN", iercod, 7L);
5870 AdvApp2Var_SysBase::mgsomsg_("MMHJCAN", 7L);
5875 //=======================================================================
5876 //function : AdvApp2Var_MathBase::mminltt_
5878 //=======================================================================
5879 int AdvApp2Var_MathBase::mminltt_(integer *ncolmx,
5885 doublereal *,//epseg,
5888 /* System generated locals */
5889 integer tabtri_dim1, tabtri_offset, i__1, i__2;
5891 /* Local variables */
5892 static logical idbg;
5893 static integer icol, ilgn, nlgn, noct, inser;
5894 static doublereal epsega;
5897 /* ***********************************************************************
5902 /* . Insert a line in a table parsed without redundance */
5906 /* TOUS,MATH_ACCES :: TABLEAU&,INSERTION,&TABLEAU */
5908 /* INPUT ARGUMENTS : */
5909 /* -------------------- */
5910 /* . NCOLMX : Number of columns in the table */
5911 /* . NLGNMX : Number of lines in the table */
5912 /* . TABTRI : Table parsed by lines without redundances */
5913 /* . NBRCOL : Number of columns used */
5914 /* . NBRLGN : Number of lines used */
5915 /* . AJOUTE : Line to be added */
5916 /* . EPSEGA : Epsilon to test the redundance */
5918 /* OUTPUT ARGUMENTS : */
5919 /* --------------------- */
5920 /* . TABTRI : Table parsed by lines without redundances */
5921 /* . NBRLGN : Number of lines used */
5922 /* . IERCOD : 0 -> No problem */
5923 /* 1 -> The table is full */
5925 /* COMMONS USED : */
5926 /* ------------------ */
5928 /* REFERENCES CALLED : */
5929 /* --------------------- */
5931 /* DESCRIPTION/NOTES/LIMITATIONS : */
5932 /* ----------------------------------- */
5933 /* . The line is inserted only if there is no line with all
5935 /* elements equl to those which are planned to be insered, to epsilon. */
5937 /* . Level of de debug = 3 */
5941 /* DECLARATIONS , CONTROL OF INPUT ARGUMENTS , INITIALIZATION */
5942 /* ***********************************************************************
5945 /* --- Parameters */
5951 /* --- Local variables */
5956 /* Parameter adjustments */
5957 tabtri_dim1 = *ncolmx;
5958 tabtri_offset = tabtri_dim1 + 1;
5959 tabtri -= tabtri_offset;
5963 ibb = AdvApp2Var_SysBase::mnfndeb_();
5966 AdvApp2Var_SysBase::mgenmsg_("MMINLTT", 7L);
5969 /* --- Control arguments */
5971 if (*nbrlgn >= *nlgnmx) {
5975 /* -------------------- */
5976 /* *** INITIALIZATION */
5977 /* -------------------- */
5981 /* ---------------------------- */
5982 /* *** SEARCH OF REDUNDANCE */
5983 /* ---------------------------- */
5986 for (ilgn = 1; ilgn <= i__1; ++ilgn) {
5987 if (tabtri[ilgn * tabtri_dim1 + 1] >= ajoute[1] - epsega) {
5988 if (tabtri[ilgn * tabtri_dim1 + 1] <= ajoute[1] + epsega) {
5990 for (icol = 1; icol <= i__2; ++icol) {
5991 if (tabtri[icol + ilgn * tabtri_dim1] < ajoute[icol] -
5992 epsega || tabtri[icol + ilgn * tabtri_dim1] >
5993 ajoute[icol] + epsega) {
6007 /* ----------------------------------- */
6008 /* *** SEARCH OF THE INSERTION POINT */
6009 /* ----------------------------------- */
6014 for (ilgn = 1; ilgn <= i__1; ++ilgn) {
6016 for (icol = 1; icol <= i__2; ++icol) {
6017 if (tabtri[icol + ilgn * tabtri_dim1] < ajoute[icol]) {
6020 if (tabtri[icol + ilgn * tabtri_dim1] > ajoute[icol]) {
6031 /* -------------- */
6033 /* -------------- */
6040 /* --- Shift lower */
6042 nlgn = *nbrlgn - inser;
6044 noct = (*ncolmx << 3) * nlgn;
6045 AdvApp2Var_SysBase::mcrfill_((integer *)&noct,
6046 (char *)&tabtri[inser * tabtri_dim1 + 1],
6047 (char *)&tabtri[(inser + 1)* tabtri_dim1 + 1]);
6052 noct = *nbrcol << 3;
6053 AdvApp2Var_SysBase::mcrfill_((integer *)&noct,
6055 (char *)&tabtri[inser * tabtri_dim1 + 1]);
6059 /* ******************************************************************** */
6060 /* OUTPUT ERROR , RETURN CALLING PROGRAM , MESSAGES */
6061 /* ******************************************************************** */
6063 /* --- The table is already full */
6072 AdvApp2Var_SysBase::maermsg_("MMINLTT", iercod, 7L);
6075 AdvApp2Var_SysBase::mgsomsg_("MMINLTT", 7L);
6080 //=======================================================================
6081 //function : AdvApp2Var_MathBase::mmjacan_
6083 //=======================================================================
6084 int AdvApp2Var_MathBase::mmjacan_(integer *ideriv,
6089 /* System generated locals */
6090 integer poljac_dim1, i__1, i__2;
6092 /* Local variables */
6093 static integer iptt, i__, j, ibb;
6094 static doublereal bid;
6096 /* ***********************************************************************
6101 /* Routine of transfer of Jacobi normalized to canonic [-1,1], */
6102 /* the tables are ranked by even, then by uneven degree. */
6106 /* LEGENDRE,JACOBI,PASSAGE. */
6108 /* INPUT ARGUMENTS : */
6109 /* ------------------ */
6110 /* IDERIV : Order of Jacobi between -1 and 2. */
6111 /* NDEG : The true degree of the polynom. */
6112 /* POLJAC : The polynom in the Jacobi base. */
6114 /* OUTPUT ARGUMENTS : */
6115 /* ------------------- */
6116 /* POLCAN : The curve expressed in the canonic base [-1,1]. */
6118 /* COMMONS USED : */
6119 /* ---------------- */
6121 /* REFERENCES CALLED : */
6122 /* ----------------------- */
6124 /* DESCRIPTION/NOTES/LIMITATIONS : */
6125 /* ----------------------------------- */
6128 /* ***********************************************************************
6131 /* Name of the routine */
6133 /* Matrices of conversion */
6136 /* ***********************************************************************
6141 /* MATRIX OF TRANSFORMATION OF LEGENDRE BASE */
6147 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
6148 /* ----------------------------------- */
6151 /* ***********************************************************************
6156 /* Legendre common / Restricted Casteljau. */
6158 /* 0:1 0 Concerns the even terms, 1 the uneven terms. */
6159 /* CANPLG : Matrix of passage to canonic from Jacobi with calculated parities */
6160 /* PLGCAN : Matrix of passage from Jacobi to canonic with calculated parities */
6163 /* ***********************************************************************
6166 /* Parameter adjustments */
6167 poljac_dim1 = *ndeg / 2 + 1;
6170 ibb = AdvApp2Var_SysBase::mnfndeb_();
6172 AdvApp2Var_SysBase::mgenmsg_("MMJACAN", 7L);
6175 /* ----------------- Expression of terms of even degree ----------------
6179 for (i__ = 0; i__ <= i__1; ++i__) {
6181 iptt = i__ * 31 - (i__ + 1) * i__ / 2 + 1;
6183 for (j = i__; j <= i__2; ++j) {
6184 bid += mmjcobi_.plgcan[iptt + j + *ideriv * 992 + 991] * poljac[
6188 polcan[i__ * 2] = bid;
6192 /* --------------- Expression of terms of uneven degree ----------------
6199 i__1 = (*ndeg - 1) / 2;
6200 for (i__ = 0; i__ <= i__1; ++i__) {
6202 iptt = i__ * 31 - (i__ + 1) * i__ / 2 + 1;
6203 i__2 = (*ndeg - 1) / 2;
6204 for (j = i__; j <= i__2; ++j) {
6205 bid += mmjcobi_.plgcan[iptt + j + ((*ideriv << 1) + 1) * 496 +
6206 991] * poljac[j + poljac_dim1];
6209 polcan[(i__ << 1) + 1] = bid;
6213 /* -------------------------------- The end -----------------------------
6218 AdvApp2Var_SysBase::mgsomsg_("MMJACAN", 7L);
6223 //=======================================================================
6224 //function : AdvApp2Var_MathBase::mmjaccv_
6226 //=======================================================================
6227 int AdvApp2Var_MathBase::mmjaccv_(integer *ncoef,
6235 /* Initialized data */
6237 static char nomprg[8+1] = "MMJACCV ";
6239 /* System generated locals */
6240 integer crvlgd_dim1, crvlgd_offset, crvcan_dim1, crvcan_offset,
6241 polaux_dim1, i__1, i__2;
6243 /* Local variables */
6244 static integer ndeg, i__, nd, ii, ibb;
6246 /* ***********************************************************************
6251 /* Passage from the normalized Jacobi base to the canonic base. */
6255 /* SMOOTHING, BASE, LEGENDRE */
6258 /* INPUT ARGUMENTS : */
6259 /* ------------------ */
6260 /* NDIM: Space Dimension. */
6261 /* NCOEF: Degree +1 of the polynom. */
6262 /* IDER: Order of Jacobi polynoms. */
6263 /* CRVLGD : Curve in the base of Jacobi. */
6265 /* OUTPUT ARGUMENTS : */
6266 /* ------------------- */
6267 /* POLAUX : Auxilliary space. */
6268 /* CRVCAN : The curve in the canonic base [-1,1] */
6270 /* COMMONS USED : */
6271 /* ---------------- */
6273 /* REFERENCES CALLED : */
6274 /* ----------------------- */
6276 /* DESCRIPTION/NOTES/LIMITATIONS : */
6277 /* ----------------------------------- */
6280 /* *********************************************************************
6283 /* Name of the routine */
6284 /* Parameter adjustments */
6285 polaux_dim1 = (*ncoef - 1) / 2 + 1;
6286 crvcan_dim1 = *ncoef - 1 + 1;
6287 crvcan_offset = crvcan_dim1;
6288 crvcan -= crvcan_offset;
6289 crvlgd_dim1 = *ncoef - 1 + 1;
6290 crvlgd_offset = crvlgd_dim1;
6291 crvlgd -= crvlgd_offset;
6295 ibb = AdvApp2Var_SysBase::mnfndeb_();
6297 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
6303 for (nd = 1; nd <= i__1; ++nd) {
6304 /* Loading of the auxilliary table. */
6307 for (i__ = 0; i__ <= i__2; ++i__) {
6308 polaux[i__] = crvlgd[ii + nd * crvlgd_dim1];
6315 i__2 = (ndeg - 1) / 2;
6316 for (i__ = 0; i__ <= i__2; ++i__) {
6317 polaux[i__ + polaux_dim1] = crvlgd[ii + nd * crvlgd_dim1];
6322 /* Call the routine of base change. */
6323 AdvApp2Var_MathBase::mmjacan_(ider, &ndeg, polaux, &crvcan[nd * crvcan_dim1]);
6332 //=======================================================================
6333 //function : mmloncv_
6335 //=======================================================================
6336 int mmloncv_(integer *ndimax,
6346 /* Initialized data */
6348 static integer kgar = 0;
6350 /* System generated locals */
6351 integer courbe_dim1, courbe_offset, i__1, i__2;
6353 /* Local variables */
6354 static doublereal tran;
6355 static integer ngaus;
6356 static doublereal c1, c2, d1, d2, wgaus[20], uroot[20], x1, x2, dd;
6357 static integer ii, jj, kk;
6358 static doublereal som;
6359 static doublereal der1, der2;
6364 /* **********************************************************************
6367 /* FUNCTION : Length of an arc of curve on a given interval */
6368 /* ---------- for a function the mathematic representation */
6369 /* which of is a multidimensional polynom. */
6370 /* The polynom is a set of polynoms the coefficients which of are ranked
6371 /* in a table with 2 indices, each line relative to 1 polynom. */
6372 /* The polynom is defined by its coefficients ordered by increasing
6373 * power of the variable. */
6374 /* All polynoms have the same number of coefficients (and the same degree). */
6376 /* KEYWORDS : LENGTH, CURVE */
6379 /* INPUT ARGUMENTS : */
6380 /* -------------------- */
6382 /* NDIMAX : Max number of lines of tables (max number of polynoms). */
6383 /* NDIMEN : Dimension of the polynom (Nomber of polynoms). */
6384 /* NCOEFF : Number of coefficients of the polynom (no limitation) */
6385 /* This is degree + 1 */
6386 /* COURBE : Coefficients of the polynom ordered by increasing power */
6387 /* Dimension to (NDIMAX,NCOEFF). */
6388 /* TDEBUT : Lower limit of integration for length calculation. */
6389 /* TFINAL : Upper limit of integration for length calculation. */
6391 /* OUTPUT ARGUMENTS : */
6392 /* --------------------- */
6393 /* XLONGC : Length of arc of curve */
6395 /* IERCOD : Error code : */
6396 /* = 0 ==> All is OK */
6397 /* = 1 ==> NDIMEN or NCOEFF negative or null */
6398 /* = 2 ==> Pb loading Legendre roots and Gauss weight */
6401 /* If error => XLONGC = 0 */
6403 /* COMMONS USED : */
6404 /* ------------------ */
6408 /* REFERENCES CALLED : */
6409 /* ---------------------- */
6411 /* MAERMSG R*8 DSQRT I*4 MIN */
6414 /* DESCRIPTION/NOTES/LIMITATIONS : */
6415 /* ----------------------------------- */
6417 /* See VGAUSS to understand well the technique. */
6418 /* Actually SQRT (dpi^2) is integrated for i=1,nbdime */
6419 /* Calculation of the derivative is included in the code to avoid an additional */
6420 /* call of the routine. */
6422 /* The integrated function is strictly increasing, it */
6423 /* is not necessary to use a high degree for the GAUSS method GAUSS. */
6425 /* The degree of LEGENDRE polynom results from the degree of the */
6426 /* polynom to be integrated. It can vary from 4 to 40 (with step of 4). */
6428 /* The precision (relative) of integration is of order 1.D-8. */
6430 /* ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
6432 /* Attention : the precision of the result is not controlled. */
6433 /* If you wish to control it, use MMCGLC1, taking into account that */
6434 /* the performance (in time) will be worse. */
6436 /* >=====================================================================
6439 /* ATTENTION : SAVE KGAR WGAUS and UROOT EVENTUALLY */
6441 /* INTEGER I1,I20 */
6442 /* PARAMETER (I1=1,I20=20) */
6444 /* Parameter adjustments */
6445 courbe_dim1 = *ndimax;
6446 courbe_offset = courbe_dim1 + 1;
6447 courbe -= courbe_offset;
6451 /* ****** General initialization ** */
6456 /* ****** Initialization of UROOT, WGAUS, NGAUS and KGAR ** */
6458 /* CALL MXVINIT(IERXV,'INTEGER',I1,KGAR,'INTEGER',I1,NGAUS */
6459 /* 1 ,'DOUBLE PRECISION',I20,UROOT,'DOUBLE PRECISION',I20,WGAUS) */
6460 /* IF (IERXV.GT.0) KGAR=0 */
6462 /* ****** Test the equity of limits ** */
6464 if (*tdebut == *tfinal) {
6469 /* ****** Test the dimension and the number of coefficients ** */
6471 if (*ndimen <= 0 || *ncoeff <= 0) {
6476 /* ****** Calculate the optimal degree ** */
6478 kk = *ncoeff / 4 + 1;
6479 kk = advapp_min(kk,10);
6481 /* ****** Return the coefficients for the integral (DEGRE=4*KK) */
6482 /* if KK <> KGAR. */
6485 mvgaus0_(&kk, uroot, wgaus, &ngaus, iercod);
6494 /* C1 => Point medium interval */
6495 /* C2 => 1/2 amplitude interval */
6497 c1 = (*tfinal + *tdebut) * .5;
6498 c2 = (*tfinal - *tdebut) * .5;
6500 /* ----------------------------------------------------------- */
6501 /* ****** Integration - Loop on GAUSS intervals ** */
6502 /* ----------------------------------------------------------- */
6507 for (jj = 1; jj <= i__1; ++jj) {
6509 /* ****** Integration taking the symmetry into account ** */
6511 tran = c2 * uroot[jj - 1];
6515 /* ****** Derivation on the dimension of the space ** */
6520 for (kk = 1; kk <= i__2; ++kk) {
6521 d1 = (*ncoeff - 1) * courbe[kk + *ncoeff * courbe_dim1];
6523 for (ii = *ncoeff - 1; ii >= 2; --ii) {
6524 dd = (ii - 1) * courbe[kk + ii * courbe_dim1];
6534 /* ****** Integration ** */
6536 som += wgaus[jj - 1] * c2 * (sqrt(der1) + sqrt(der2));
6538 /* ****** End of loop on GAUSS intervals ** */
6543 /* ****** Work ended ** */
6547 /* ****** It is forced IERCOD = 0 ** */
6551 /* ****** Final processing ** */
6555 /* ****** Save UROOT, WGAUS, NGAUS and KGAR ** */
6557 /* CALL MXVSAVE(IERXV,'INTEGER',I1,KGAR,'INTEGER',I1,NGAUS */
6558 /* 1 ,'DOUBLE PRECISION',I20,UROOT,'DOUBLE PRECISION',I20,WGAUS) */
6559 /* IF (IERXV.GT.0) KGAR=0 */
6561 /* ****** End of sub-program ** */
6564 AdvApp2Var_SysBase::maermsg_("MMLONCV", iercod, 7L);
6569 //=======================================================================
6570 //function : AdvApp2Var_MathBase::mmpobas_
6572 //=======================================================================
6573 int AdvApp2Var_MathBase::mmpobas_(doublereal *tparam,
6581 static integer c__2 = 2;
6582 static integer c__1 = 1;
6585 /* Initialized data */
6587 static doublereal moin11[2] = { -1.,1. };
6589 /* System generated locals */
6590 integer valbas_dim1, i__1;
6592 /* Local variables */
6593 static doublereal vjac[80], herm[24];
6594 static integer iord[2];
6595 static doublereal wval[4];
6596 static integer nwcof, iunit;
6597 static doublereal wpoly[7];
6598 static integer ii, jj, iorjac;
6599 static doublereal hermit[36] /* was [6][3][2] */;
6600 static integer kk1, kk2, kk3;
6601 static integer khe, ier;
6604 /* ***********************************************************************
6609 /* Position on the polynoms of base hermit-Jacobi */
6610 /* and their succesive derivatives */
6614 /* PUBLIC, POSITION, HERMIT, JACOBI */
6616 /* INPUT ARGUMENTS : */
6617 /* -------------------- */
6618 /* TPARAM : Parameter for which the position is found. */
6619 /* IORDRE : Orderof hermit-Jacobi (-1,0,1, ou 2) */
6620 /* NCOEFF : Number of coefficients of polynoms (Nb of value to calculate) */
6621 /* NDERIV : Number of derivative to calculate (0<= N <=3) */
6622 /* 0 -> Position simple on base functions */
6623 /* N -> Position on base functions and derivative */
6624 /* of order 1 to N */
6626 /* OUTPUT ARGUMENTS : */
6627 /* --------------------- */
6628 /* VALBAS (NCOEFF, 0:NDERIV) : calculated value */
6630 /* d vj(t) = VALBAS(J, I) */
6634 /* IERCOD : Error code */
6636 /* 1 : Incoherence of input arguments */
6638 /* COMMONS USED : */
6639 /* -------------- */
6642 /* REFERENCES CALLED : */
6643 /* ------------------- */
6646 /* DESCRIPTION/NOTES/LIMITATIONS : */
6647 /* ----------------------------------- */
6650 /* ***********************************************************************
6653 /* ***********************************************************************
6658 /* Parameter adjustments */
6659 valbas_dim1 = *ncoeff;
6664 /* ***********************************************************************
6666 /* INITIALIZATIONS */
6667 /* ***********************************************************************
6672 /* ***********************************************************************
6675 /* ***********************************************************************
6690 iorjac = (*iordre + 1) << 1;
6692 /* (1) Generic Calculations .... */
6694 /* (1.a) Calculation of hermit polynoms */
6697 mmherm1_(moin11, &c__2, iord, hermit, &ier);
6703 /* (1.b) Evaluation of hermit polynoms */
6706 iunit = *nderiv + 1;
6707 khe = (*iordre + 1) * iunit;
6712 for (ii = 0; ii <= i__1; ++ii) {
6713 mmdrvcb_(nderiv, &c__1, &iorjac, &hermit[(ii + 3) * 6 - 18],
6714 tparam, &herm[jj - 1], &ier);
6719 mmdrvcb_(nderiv, &c__1, &iorjac, &hermit[(ii + 6) * 6 - 18],
6720 tparam, &herm[jj + khe - 1], &ier);
6730 for (ii = 0; ii <= i__1; ++ii) {
6731 AdvApp2Var_MathBase::mmpocrb_(&c__1, &iorjac, &hermit[(ii + 3) * 6 - 18], &c__1,
6732 tparam, &herm[jj - 1]);
6734 AdvApp2Var_MathBase::mmpocrb_(&c__1, &iorjac, &hermit[(ii + 6) * 6 - 18], &c__1,
6735 tparam, &herm[jj + khe - 1]);
6740 /* (1.c) Evaluation of Jacobi polynoms */
6742 ii = *ncoeff - iorjac;
6744 mmpojac_(tparam, &iorjac, &ii, nderiv, vjac, &ier);
6749 /* (1.d) Evaluation of W(t) */
6753 nwcof = advapp_max(i__1,1);
6754 AdvApp2Var_SysBase::mvriraz_((integer *)&nwcof,
6761 } else if (*iordre == 1) {
6764 } else if (*iordre == 0) {
6768 mmdrvcb_(nderiv, &c__1, &nwcof, wpoly, tparam, wval, &ier);
6773 kk1 = *ncoeff - iorjac;
6777 /* (2) Evaluation of order 0 */
6781 for (ii = 1; ii <= i__1; ++ii) {
6782 valbas[ii] = herm[jj - 1];
6787 for (ii = 1; ii <= i__1; ++ii) {
6788 valbas[ii + iorjac] = wval[0] * vjac[ii - 1];
6791 /* (3) Evaluation of order 1 */
6796 for (ii = 1; ii <= i__1; ++ii) {
6797 valbas[ii + valbas_dim1] = herm[jj - 1];
6803 for (ii = 1; ii <= i__1; ++ii) {
6804 valbas[ii + iorjac + valbas_dim1] = wval[0] * vjac[ii + kk1 - 1]
6805 + wval[1] * vjac[ii - 1];
6809 /* (4) Evaluation of order 2 */
6814 for (ii = 1; ii <= i__1; ++ii) {
6815 valbas[ii + (valbas_dim1 << 1)] = herm[jj - 1];
6820 for (ii = 1; ii <= i__1; ++ii) {
6821 valbas[ii + iorjac + (valbas_dim1 << 1)] = wval[0] * vjac[ii +
6822 kk2 - 1] + wval[1] * 2 * vjac[ii + kk1 - 1] + wval[2] *
6827 /* (5) Evaluation of order 3 */
6832 for (ii = 1; ii <= i__1; ++ii) {
6833 valbas[ii + valbas_dim1 * 3] = herm[jj - 1];
6838 for (ii = 1; ii <= i__1; ++ii) {
6839 valbas[ii + iorjac + valbas_dim1 * 3] = wval[0] * vjac[ii + kk3 -
6840 1] + wval[1] * 3 * vjac[ii + kk2 - 1] + wval[2] * 3 *
6841 vjac[ii + kk1 - 1] + wval[3] * vjac[ii - 1];
6847 /* ***********************************************************************
6849 /* ERROR PROCESSING */
6850 /* ***********************************************************************
6860 /* ***********************************************************************
6862 /* RETURN CALLING PROGRAM */
6863 /* ***********************************************************************
6869 AdvApp2Var_SysBase::maermsg_("MMPOBAS", iercod, 7L);
6874 //=======================================================================
6875 //function : AdvApp2Var_MathBase::mmpocrb_
6877 //=======================================================================
6878 int AdvApp2Var_MathBase::mmpocrb_(integer *ndimax,
6886 /* System generated locals */
6887 integer courbe_dim1, courbe_offset, i__1, i__2;
6889 /* Local variables */
6890 static integer ncof2;
6891 static integer isize, nd, kcf, ncf;
6894 /* ***********************************************************************
6899 /* CALCULATE THE COORDINATES OF A POINT OF A CURVE OF GIVEN PARAMETER */
6900 /* TPARAM ( IN 2D, 3D OR MORE) */
6904 /* TOUS , MATH_ACCES :: COURBE&,PARAMETRE& , POSITIONNEMENT , &POINT
6907 /* INPUT ARGUMENTS : */
6908 /* ------------------ */
6909 /* NDIMAX : format / dimension of the curve */
6910 /* NCOEFF : Nb of coefficients of the curve */
6911 /* COURBE : Matrix of coefficients of the curve */
6912 /* NDIM : Dimension useful of the workspace */
6913 /* TPARAM : Value of the parameter where the point is calculated */
6915 /* OUTPUT ARGUMENTS : */
6916 /* ------------------- */
6917 /* PNTCRB : Coordinates of the calculated point */
6919 /* COMMONS USED : */
6920 /* ---------------- */
6924 /* REFERENCES CALLED : */
6925 /* ---------------------- */
6927 /* MIRAZ MVPSCR2 MVPSCR3 */
6929 /* DESCRIPTION/NOTES/LIMITATIONS : */
6930 /* ----------------------------------- */
6933 /* ***********************************************************************
6937 /* ***********************************************************************
6940 /* Parameter adjustments */
6941 courbe_dim1 = *ndimax;
6942 courbe_offset = courbe_dim1 + 1;
6943 courbe -= courbe_offset;
6948 AdvApp2Var_SysBase::miraz_((integer *)&isize,
6949 (char *)&pntcrb[1]);
6955 /* optimal processing 3d */
6957 if (*ndim == 3 && *ndimax == 3) {
6958 mvpscr3_(ncoeff, &courbe[courbe_offset], tparam, &pntcrb[1]);
6960 /* optimal processing 2d */
6962 } else if (*ndim == 2 && *ndimax == 2) {
6963 mvpscr2_(ncoeff, &courbe[courbe_offset], tparam, &pntcrb[1]);
6965 /* Any dimension - scheme of HORNER */
6967 } else if (*tparam == 0.) {
6969 for (nd = 1; nd <= i__1; ++nd) {
6970 pntcrb[nd] = courbe[nd + courbe_dim1];
6973 } else if (*tparam == 1.) {
6975 for (ncf = 1; ncf <= i__1; ++ncf) {
6977 for (nd = 1; nd <= i__2; ++nd) {
6978 pntcrb[nd] += courbe[nd + ncf * courbe_dim1];
6984 ncof2 = *ncoeff + 2;
6986 for (nd = 1; nd <= i__1; ++nd) {
6988 for (ncf = 2; ncf <= i__2; ++ncf) {
6990 pntcrb[nd] = (pntcrb[nd] + courbe[nd + kcf * courbe_dim1]) * *
6994 pntcrb[nd] += courbe[nd + courbe_dim1];
7003 //=======================================================================
7004 //function : AdvApp2Var_MathBase::mmmpocur_
7006 //=======================================================================
7007 int AdvApp2Var_MathBase::mmmpocur_(integer *ncofmx,
7015 /* System generated locals */
7016 integer courbe_dim1, courbe_offset, i__1;
7018 /* Local variables */
7019 static integer i__, nd;
7020 static doublereal fu;
7023 /* ***********************************************************************
7028 /* Position of a point on curve (ncofmx,ndim). */
7032 /* TOUS , AB_SPECIFI :: COURBE&,POLYNOME&,POSITIONNEMENT,&POINT */
7034 /* INPUT ARGUMENTS : */
7035 /* ------------------ */
7036 /* NCOFMX: Format / degree of the CURVE. */
7037 /* NDIM : Dimension of the space. */
7038 /* NDEG : Degree of the polynom. */
7039 /* COURBE: Coefficients of the curve. */
7040 /* TPARAM: Parameter on the curve */
7042 /* OUTPUT ARGUMENTS : */
7043 /* ------------------- */
7044 /* TABVAL(NDIM): The resulting point (or table of values) */
7046 /* COMMONS USED : */
7047 /* ---------------- */
7049 /* REFERENCES CALLED : */
7050 /* ----------------------- */
7052 /* DESCRIPTION/NOTES/LIMITATIONS : */
7053 /* ----------------------------------- */
7056 /* ***********************************************************************
7059 /* Parameter adjustments */
7061 courbe_dim1 = *ncofmx;
7062 courbe_offset = courbe_dim1 + 1;
7063 courbe -= courbe_offset;
7068 for (nd = 1; nd <= i__1; ++nd) {
7074 for (nd = 1; nd <= i__1; ++nd) {
7075 fu = courbe[*ndeg + nd * courbe_dim1];
7076 for (i__ = *ndeg - 1; i__ >= 1; --i__) {
7077 fu = fu * *tparam + courbe[i__ + nd * courbe_dim1];
7087 //=======================================================================
7088 //function : mmpojac_
7090 //=======================================================================
7091 int mmpojac_(doublereal *tparam,
7099 static integer c__2 = 2;
7101 /* Initialized data */
7103 static integer nbcof = -1;
7105 /* System generated locals */
7106 integer valjac_dim1, i__1, i__2;
7108 /* Local variables */
7109 static doublereal cofa, cofb, denom, tnorm[100];
7110 static integer ii, jj, kk1, kk2;
7111 static doublereal aux1, aux2;
7114 /* ***********************************************************************
7119 /* Positioning on Jacobi polynoms and their derivatives */
7120 /* successive by a recurrent algorithm */
7124 /* RESERVE, POSITIONING, JACOBI */
7126 /* INPUT ARGUMENTS : */
7127 /* -------------------- */
7128 /* TPARAM : Parameter for which positioning is done. */
7129 /* IORDRE : Order of hermit-?? (-1,0,1, or 2) */
7130 /* NCOEFF : Number of coeeficients of polynoms (Nb of value to */
7132 /* NDERIV : Number of derivative to calculate (0<= N <=3) */
7133 /* 0 -> Position simple on jacobi functions */
7134 /* N -> Position on jacobi functions and their */
7135 /* derivatives of order 1 to N. */
7137 /* OUTPUT ARGUMENTS : */
7138 /* --------------------- */
7139 /* VALJAC (NCOEFF, 0:NDERIV) : the calculated values */
7141 /* d vj(t) = VALJAC(J, I) */
7145 /* IERCOD : Error Code */
7147 /* 1 : Incoherence of input arguments */
7149 /* COMMONS USED : */
7150 /* ------------------ */
7153 /* REFERENCES CALLED : */
7154 /* --------------------- */
7157 /* DESCRIPTION/NOTES/LIMITATIONS : */
7158 /* ----------------------------------- */
7161 /* ***********************************************************************
7164 /* ***********************************************************************
7168 /* static varaibles */
7172 /* Parameter adjustments */
7173 valjac_dim1 = *ncoeff;
7178 /* ***********************************************************************
7180 /* INITIALISATIONS */
7181 /* ***********************************************************************
7186 /* ***********************************************************************
7189 /* ***********************************************************************
7195 if (*ncoeff > 100) {
7199 /* --- Calculation of norms */
7201 /* IF (NCOEFF.GT.NBCOF) THEN */
7203 for (ii = 1; ii <= i__1; ++ii) {
7207 for (jj = 1; jj <= i__2; ++jj) {
7208 aux2 = aux2 * (doublereal) (kk1 + *iordre + jj) / (doublereal) (
7211 i__2 = (*iordre << 1) + 1;
7212 tnorm[ii - 1] = sqrt(aux2 * (kk1 * 2. + (*iordre << 1) + 1) / pow__ii(&
7220 /* --- Trivial Positions ----- */
7223 aux1 = (doublereal) (*iordre + 1);
7224 valjac[2] = aux1 * *tparam;
7227 valjac[valjac_dim1 + 1] = 0.;
7228 valjac[valjac_dim1 + 2] = aux1;
7231 valjac[(valjac_dim1 << 1) + 1] = 0.;
7232 valjac[(valjac_dim1 << 1) + 2] = 0.;
7235 valjac[valjac_dim1 * 3 + 1] = 0.;
7236 valjac[valjac_dim1 * 3 + 2] = 0.;
7241 /* --- Positioning by recurrence */
7244 for (ii = 3; ii <= i__1; ++ii) {
7248 aux1 = (doublereal) (*iordre + kk2);
7250 cofa = aux2 * (aux2 + 1) * (aux2 + 2);
7251 cofb = (aux2 + 2) * -2. * aux1 * aux1;
7252 denom = kk1 * 2. * (kk2 + (*iordre << 1) + 1) * aux2;
7256 valjac[ii] = (cofa * *tparam * valjac[kk1] + cofb * valjac[kk2]) *
7260 valjac[ii + valjac_dim1] = (cofa * *tparam * valjac[kk1 +
7261 valjac_dim1] + cofa * valjac[kk1] + cofb * valjac[kk2 +
7262 valjac_dim1]) * denom;
7265 valjac[ii + (valjac_dim1 << 1)] = (cofa * *tparam * valjac[
7266 kk1 + (valjac_dim1 << 1)] + cofa * 2 * valjac[kk1 +
7267 valjac_dim1] + cofb * valjac[kk2 + (valjac_dim1 << 1)]
7272 valjac[ii + valjac_dim1 * 3] = (cofa * *tparam * valjac[kk1 +
7273 valjac_dim1 * 3] + cofa * 3 * valjac[kk1 + (
7274 valjac_dim1 << 1)] + cofb * valjac[kk2 + valjac_dim1 *
7280 /* ---> Normalization */
7283 for (ii = 1; ii <= i__1; ++ii) {
7285 for (jj = 0; jj <= i__2; ++jj) {
7286 valjac[ii + jj * valjac_dim1] = tnorm[ii - 1] * valjac[ii + jj *
7293 /* ***********************************************************************
7295 /* PROCESSING OF ERRORS */
7296 /* ***********************************************************************
7304 /* ***********************************************************************
7306 /* RETURN CALLING PROGRAM */
7307 /* ***********************************************************************
7313 AdvApp2Var_SysBase::maermsg_("MMPOJAC", iercod, 7L);
7318 //=======================================================================
7319 //function : AdvApp2Var_MathBase::mmposui_
7321 //=======================================================================
7322 int AdvApp2Var_MathBase::mmposui_(integer *dimmat,
7329 /* System generated locals */
7332 /* Local variables */
7333 static logical ldbg;
7334 static integer imin, jmin, i__, j, k;
7335 static logical trouve;
7337 /* ***********************************************************************
7342 /* FILL THE TABLE OF POSITIONING POSUIV WHICH ALLOWS TO */
7343 /* PARSE BY COLUMN THE INFERIOR TRIANGULAR PART OF THE */
7344 /* MATRIX IN FORM OF PROFILE */
7349 /* RESERVE, MATRIX, PROFILE */
7351 /* INPUT ARGUMENTS : */
7352 /* -------------------- */
7354 /* NISTOC: NUMBER OF COEFFICIENTS IN THE PROFILE */
7355 /* DIMMAT: NUMBER OF LINE OF THE SYMMETRIC SQUARE MATRIX */
7356 /* APOSIT: TABLE OF POSITIONING OF STORAGE TERMS */
7357 /* APOSIT(1,I) CONTAINS THE NUMBER OF TERMES-1 ON LINE
7358 /* I IN THE PROFILE OF THE MATRIX */
7359 /* APOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF DIAGONAL TERM
7363 /* OUTPUT ARGUMENTS : */
7364 /* --------------------- */
7365 /* POSUIV: POSUIV(K) (WHERE K IS THE INDEX OF STORAGE OF MAT(I,J)) */
7366 /* CONTAINS THE SMALLEST NUMBER IMIN>I OF THE LINE THAT */
7367 /* POSSESSES A TERM MAT(IMIN,J) THAT IS IN THE PROFILE. */
7368 /* IF THERE IS NO TERM MAT(IMIN,J) IN THE PROFILE THEN POSUIV(K)=-1 */
7371 /* COMMONS USED : */
7372 /* ------------------ */
7375 /* REFERENCES CALLED : */
7376 /* --------------------- */
7379 /* DESCRIPTION/NOTES/LIMITATIONS : */
7380 /* ----------------------------------- */
7383 /* ***********************************************************************
7386 /* ***********************************************************************
7391 /* ***********************************************************************
7393 /* INITIALIZATIONS */
7394 /* ***********************************************************************
7397 /* Parameter adjustments */
7402 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
7404 AdvApp2Var_SysBase::mgenmsg_("MMPOSUI", 7L);
7409 /* ***********************************************************************
7412 /* ***********************************************************************
7418 for (i__ = 1; i__ <= i__1; ++i__) {
7419 jmin = i__ - aposit[(i__ << 1) + 1];
7421 for (j = jmin; j <= i__2; ++j) {
7424 while(! trouve && imin <= *dimmat) {
7425 if (imin - aposit[(imin << 1) + 1] <= j) {
7431 k = aposit[(i__ << 1) + 2] - i__ + j;
7446 /* ***********************************************************************
7448 /* ERROR PROCESSING */
7449 /* ***********************************************************************
7455 /* ***********************************************************************
7457 /* RETURN CALLING PROGRAM */
7458 /* ***********************************************************************
7463 /* ___ DESALLOCATION, ... */
7465 AdvApp2Var_SysBase::maermsg_("MMPOSUI", iercod, 7L);
7467 AdvApp2Var_SysBase::mgsomsg_("MMPOSUI", 7L);
7472 //=======================================================================
7473 //function : AdvApp2Var_MathBase::mmresol_
7475 //=======================================================================
7476 int AdvApp2Var_MathBase::mmresol_(integer *hdimen,
7494 static integer c__100 = 100;
7496 /* System generated locals */
7499 /* Local variables */
7500 static logical ldbg;
7501 static doublereal mcho[100];
7502 static integer jmin, jmax, i__, j, k, l;
7503 static long int iofv1, iofv2, iofv3, iofv4;
7504 static doublereal v1[100], v2[100], v3[100], v4[100];
7505 static integer deblig, dimhch;
7506 static doublereal hchole[100];
7507 static long int iofmch, iofmam, iofhch;
7508 static doublereal matsym[100];
7514 /* ***********************************************************************
7519 /* SOLUTION OF THE SYSTEM */
7526 /* RESERVE, SOLUTION, SYSTEM, LAGRANGIAN */
7528 /* INPUT ARGUMENTS : */
7529 /* -------------------- */
7530 /* HDIMEN: NOMBER OF LINE (OR COLUMN) OF THE HESSIAN MATRIX */
7531 /* GDIMEN: NOMBER OF LINE OF THE MATRIX OF CONSTRAINTS */
7532 /* HNSTOC: NOMBErS OF TERMS IN THE PROFILE OF HESSIAN MATRIX
7534 /* GNSTOC: NOMBERS OF TERMS IN THE PROFILE OF THE MATRIX OF CONSTRAINTS */
7535 /* MNSTOC: NOMBERS OF TERMS IN THE PROFILE OF THE MATRIX M= G H t(G) */
7536 /* where H IS THE HESSIAN MATRIX AND G IS THE MATRIX OF CONSTRAINTS */
7537 /* MATSYH: TRIANGULAR INFERIOR PART OF THE HESSIAN MATRIX
7538 /* IN FORM OF PROFILE */
7539 /* MATSYG: MATRIX OF CONSTRAINTS IN FORM OF PROFILE */
7540 /* VECSYH: VECTOR OF THE SECOND MEMBER ASSOCIATED TO MATSYH */
7541 /* VECSYG: VECTOR OF THE SECOND MEMBER ASSOCIATED TO MATSYG */
7542 /* HPOSIT: TABLE OF POSITIONING OF THE HESSIAN MATRIX */
7543 /* HPOSIT(1,I) CONTAINS THE NUMBER OF TERMS -1 */
7544 /* WHICH ARE IN THE PROFILE AT LINE I */
7545 /* HPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF TERM */
7546 /* DIAGONAL OF THE MATRIX AT LINE I */
7547 /* HPOSUI: TABLE ALLOWING TO PARSE THE HESSIAN MATRIX BY COLUMN */
7548 /* IN FORM OF PROFILE */
7549 /* HPOSUI(K) CONTAINS THE NUMBER OF LINE IMIN FOLLOWING THE CURRENT LINE*/
7550 /* I WHERE H(I,J)=MATSYH(K) AS IT EXISTS IN THE */
7551 /* SAME COLUMN J A TERM IN THE PROFILE OF LINE IMIN */
7552 /* IF SUCH TERM DOES NOT EXIST IMIN=-1 */
7553 /* GPOSIT: TABLE OF POSITIONING OF THE MATRIX OF CONSTRAINTS */
7554 /* GPOSIT(1,I) CONTAINS THE NUMBER OF TERMS OF LINE I */
7555 /* WHICH ARE IN THE PROFILE */
7556 /* GPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF THE LAST TERM
7557 /* OF LINE I WHICH IS IN THE PROFILE */
7558 /* GPOSIT(3,I) CONTAINS THE NUMBER OF COLUMN CORRESPONDING */
7559 /* TO THE FIRST TERM OF LINE I WHICH IS IN THE PROFILE */
7560 /* MMPOSUI, MPOSIT: SAME STRUCTURE AS HPOSUI, BUT FOR MATRIX
7564 /* OUTPUT ARGUMENTS : */
7565 /* --------------------- */
7566 /* VECSOL: VECTOR SOLUTION V OF THE SYSTEM */
7567 /* IERCOD: ERROR CODE */
7569 /* COMMONS USED : */
7570 /* ------------------ */
7573 /* REFERENCES CALLED : */
7574 /* --------------------- */
7577 /* DESCRIPTION/NOTES/LIMITATIONS : */
7578 /* ----------------------------------- */
7580 /* ***********************************************************************
7583 /* ***********************************************************************
7586 /* ***********************************************************************
7588 /* INITIALISATIONS */
7589 /* ***********************************************************************
7592 /* Parameter adjustments */
7605 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
7607 AdvApp2Var_SysBase::mgenmsg_("MMRESOL", 7L);
7618 /* ***********************************************************************
7621 /* ***********************************************************************
7624 /* Dynamic allocation */
7626 AdvApp2Var_SysBase::macrar8_(hdimen, &c__100, v1, &iofv1, &ier);
7630 dimhch = hposit[(*hdimen << 1) + 2];
7631 AdvApp2Var_SysBase::macrar8_(&dimhch, &c__100, hchole, &iofhch, &ier);
7636 /* solution of system 1 H V1 = b */
7637 /* where H=MATSYH and b=VECSYH */
7639 mmchole_(hnstoc, hdimen, &matsyh[1], &hposit[3], &hposui[1], &hchole[
7644 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1], &vecsyh[
7645 1], &v1[iofv1], &ier);
7650 /* Case when there are constraints */
7654 /* Calculate the vector of the second member V2=G H(-1) b -c = G v1-c */
7655 /* of system of unknown Lagrangian vector MULTIP */
7656 /* where G=MATSYG */
7659 AdvApp2Var_SysBase::macrar8_(gdimen, &c__100, v2, &iofv2, &ier);
7663 AdvApp2Var_SysBase::macrar8_(hdimen, &c__100, v3, &iofv3, &ier);
7667 AdvApp2Var_SysBase::macrar8_(gdimen, &c__100, v4, &iofv4, &ier);
7671 AdvApp2Var_SysBase::macrar8_(mnstoc, &c__100, matsym, &iofmam, &ier);
7677 mmatvec_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v1[iofv1], &
7678 deblig, &v2[iofv2], &ier);
7683 for (i__ = 1; i__ <= i__1; ++i__) {
7684 v2[i__ + iofv2 - 1] -= vecsyg[i__];
7687 /* Calculate the matrix M= G H(-1) t(G) */
7688 /* RESOL DU SYST 2 : H qi = gi */
7689 /* where is a vector column of t(G) */
7691 /* then calculate G qi */
7692 /* then construct M in form of profile */
7697 for (i__ = 1; i__ <= i__1; ++i__) {
7698 AdvApp2Var_SysBase::mvriraz_((integer *)hdimen, (char *)&v1[iofv1]);
7699 AdvApp2Var_SysBase::mvriraz_((integer *)hdimen, (char *)&v3[iofv3]);
7700 AdvApp2Var_SysBase::mvriraz_((integer *)gdimen, (char *)&v4[iofv4]);
7701 jmin = gposit[i__ * 3 + 3];
7702 jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
7703 aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
7705 for (j = jmin; j <= i__2; ++j) {
7707 v1[j + iofv1 - 1] = matsyg[k];
7709 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1],
7710 &v1[iofv1], &v3[iofv3], &ier);
7716 mmatvec_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v3[
7717 iofv3], &deblig, &v4[iofv4], &ier);
7722 k = mposit[(i__ << 1) + 2];
7723 matsym[k + iofmam - 1] = v4[i__ + iofv4 - 1];
7724 while(mmposui[k] > 0) {
7726 k = mposit[(l << 1) + 2] - l + i__;
7727 matsym[k + iofmam - 1] = v4[l + iofv4 - 1];
7732 /* SOLVE SYST 3 M L = V2 */
7736 AdvApp2Var_SysBase::mvriraz_((integer *)gdimen, (char *)&v4[iofv4]);
7737 AdvApp2Var_SysBase::macrar8_(mnstoc, &c__100, mcho, &iofmch, &ier);
7741 mmchole_(mnstoc, gdimen, &matsym[iofmam], &mposit[3], &mmposui[1], &
7742 mcho[iofmch], &ier);
7746 mmrslss_(mnstoc, gdimen, &mcho[iofmch], &mposit[3], &mmposui[1], &v2[
7747 iofv2], &v4[iofv4], &ier);
7753 /* CALCULATE THE VECTOR OF THE SECOND MEMBER OF THE SYSTEM Hx = b - t(G) L
7757 AdvApp2Var_SysBase::mvriraz_((integer *)hdimen, (char *)&v1[iofv1]);
7758 mmtmave_(gdimen, hdimen, &gposit[4], gnstoc, &matsyg[1], &v4[iofv4], &
7764 for (i__ = 1; i__ <= i__1; ++i__) {
7765 v1[i__ + iofv1 - 1] = vecsyh[i__] - v1[i__ + iofv1 - 1];
7768 /* RESOL SYST 4 Hx = b - t(G) L */
7771 mmrslss_(hnstoc, hdimen, &hchole[iofhch], &hposit[3], &hposui[1], &v1[
7772 iofv1], &vecsol[1], &ier);
7778 for (i__ = 1; i__ <= i__1; ++i__) {
7779 vecsol[i__] = v1[i__ + iofv1 - 1];
7785 /* ***********************************************************************
7787 /* PROCESSING OF ERRORS */
7788 /* ***********************************************************************
7797 AdvApp2Var_SysBase::mswrdbg_("MMRESOL : PROBLEM WITH DIMMAT", 30L);
7800 /* ***********************************************************************
7802 /* RETURN CALLING PROGRAM */
7803 /* ***********************************************************************
7808 /* ___ DESALLOCATION, ... */
7809 AdvApp2Var_SysBase::macrdr8_(hdimen, &c__100, v1, &iofv1, &ier);
7810 if (*iercod == 0 && ier > 0) {
7813 AdvApp2Var_SysBase::macrdr8_(&dimhch, &c__100, hchole, &iofhch, &ier);
7814 if (*iercod == 0 && ier > 0) {
7817 AdvApp2Var_SysBase::macrdr8_(gdimen, &c__100, v2, &iofv2, &ier);
7818 if (*iercod == 0 && ier > 0) {
7821 AdvApp2Var_SysBase::macrdr8_(hdimen, &c__100, v3, &iofv3, &ier);
7822 if (*iercod == 0 && ier > 0) {
7825 AdvApp2Var_SysBase::macrdr8_(gdimen, &c__100, v4, &iofv4, &ier);
7826 if (*iercod == 0 && ier > 0) {
7829 AdvApp2Var_SysBase::macrdr8_(mnstoc, &c__100, matsym, &iofmam, &ier);
7830 if (*iercod == 0 && ier > 0) {
7833 AdvApp2Var_SysBase::macrdr8_(mnstoc, &c__100, mcho, &iofmch, &ier);
7834 if (*iercod == 0 && ier > 0) {
7838 AdvApp2Var_SysBase::maermsg_("MMRESOL", iercod, 7L);
7840 AdvApp2Var_SysBase::mgsomsg_("MMRESOL", 7L);
7845 //=======================================================================
7846 //function : mmrslss_
7848 //=======================================================================
7849 int mmrslss_(integer *,//mxcoef,
7854 doublereal *mscnmbr,
7858 /* System generated locals */
7861 /* Local variables */
7862 static logical ldbg;
7863 static integer i__, j;
7864 static doublereal somme;
7865 static integer pointe, ptcour;
7867 /* ***********************************************************************
7872 /* Solves linear system SS x = b where S is a */
7873 /* triangular lower matrix given in form of profile */
7877 /* RESERVE, MATRICE_PROFILE, RESOLUTION, CHOLESKI */
7879 /* INPUT ARGUMENTS : */
7880 /* -------------------- */
7881 /* MXCOEF : Maximum number of non-null coefficient in the matrix */
7882 /* DIMENS : Dimension of the matrix */
7883 /* SMATRI(MXCOEF) : Values of coefficients of the matrix */
7884 /* SPOSIT(2,DIMENS): */
7885 /* SPOSIT(1,*) : Distance diagonal-extremity of the line */
7886 /* SPOSIT(2,*) : Position of diagonal terms in AMATRI */
7887 /* POSUIV(MXCOEF): first line inferior not out of profile */
7888 /* MSCNMBR(DIMENS): Vector second member of the equation */
7890 /* OUTPUT ARGUMENTS : */
7891 /* --------------------- */
7892 /* SOLUTI(NDIMEN) : Result vector */
7893 /* IERCOD : Error code 0 : ok */
7895 /* COMMONS USED : */
7896 /* ------------------ */
7899 /* REFERENCES CALLED : */
7900 /* --------------------- */
7903 /* DESCRIPTION/NOTES/LIMITATIONS : */
7904 /* ----------------------------------- */
7906 /* SS is the decomposition of choleski of a symmetric matrix */
7907 /* defined postive, that can result from routine MMCHOLE. */
7909 /* For a full matrix it is possible to use MRSLMSC */
7911 /* LEVEL OF DEBUG = 4 */
7913 /* ***********************************************************************
7916 /* ***********************************************************************
7921 /* ***********************************************************************
7923 /* INITIALISATIONS */
7924 /* ***********************************************************************
7927 /* Parameter adjustments */
7935 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
7937 AdvApp2Var_SysBase::mgenmsg_("MMRSLSS", 7L);
7941 /* ***********************************************************************
7944 /* ***********************************************************************
7947 /* ----- Solution of Sw = b */
7950 for (i__ = 1; i__ <= i__1; ++i__) {
7952 pointe = sposit[(i__ << 1) + 2];
7955 for (j = i__ - sposit[(i__ << 1) + 1]; j <= i__2; ++j) {
7956 somme += smatri[pointe - (i__ - j)] * soluti[j];
7959 soluti[i__] = (mscnmbr[i__] - somme) / smatri[pointe];
7962 /* ----- Solution of S u = w */
7964 for (i__ = *dimens; i__ >= 1; --i__) {
7966 pointe = sposit[(i__ << 1) + 2];
7970 ptcour = sposit[(j << 1) + 2] - (j - i__);
7971 somme += smatri[ptcour] * soluti[j];
7975 soluti[i__] = (soluti[i__] - somme) / smatri[pointe];
7980 /* ***********************************************************************
7982 /* ERROR PROCESSING */
7983 /* ***********************************************************************
7987 /* ***********************************************************************
7989 /* RETURN PROGRAM CALLING */
7990 /* ***********************************************************************
7995 AdvApp2Var_SysBase::maermsg_("MMRSLSS", iercod, 7L);
7997 AdvApp2Var_SysBase::mgsomsg_("MMRSLSS", 7L);
8002 //=======================================================================
8003 //function : mmrslw_
8005 //=======================================================================
8006 int mmrslw_(integer *normax,
8014 /* System generated locals */
8015 integer abmatr_dim1, abmatr_offset, xmatri_dim1, xmatri_offset, i__1,
8019 /* Local variables */
8020 static integer kpiv;
8021 static doublereal pivot;
8022 static integer ii, jj, kk;
8023 static doublereal akj;
8026 /* **********************************************************************
8031 /* Solution of a linear system A.x = B of N equations to N */
8032 /* unknown by Gauss method (partial pivot) or : */
8033 /* A is matrix NORDRE * NORDRE, */
8034 /* B is matrix NORDRE (lines) * NDIMEN (columns), */
8035 /* x is matrix NORDRE (lines) * NDIMEN (columns). */
8036 /* In this program, A and B are stored in matrix ABMATR */
8037 /* the lines and columns which of were inverted. ABMATR(k,j) is */
8038 /* term A(j,k) if k <= NORDRE, B(j,k-NORDRE) otherwise (see example). */
8042 /* TOUS, MATH_ACCES::EQUATION&, MATRICE&, RESOLUTION, GAUSS, &SOLUTION */
8044 /* INPUT ARGUMENTS : */
8045 /* ------------------ */
8046 /* NORMAX : Max size of the first index of XMATRI. This argument */
8047 /* serves only for the declaration of dimension of XMATRI and should be */
8048 /* above or equal to NORDRE. */
8049 /* NORDRE : Order of the matrix i.e. number of equations and */
8050 /* unknown quantities of the linear system to be solved. */
8051 /* NDIMEN : Number of the second member. */
8052 /* EPSPIV : Minimal value of a pivot. If during the calculation */
8053 /* the absolute value of the pivot is below EPSPIV, the */
8054 /* system of equations is declared singular. EPSPIV should */
8055 /* be a "small" real. */
8057 /* ABMATR(NORDRE+NDIMEN,NORDRE) : Auxiliary matrix containing */
8058 /* matrix A and matrix B. */
8060 /* OUTPUT ARGUMENTS : */
8061 /* ------------------- */
8062 /* XMATRI : Matrix containing NORDRE*NDIMEN solutions. */
8063 /* IERCOD=0 shows that all solutions are calculated. */
8064 /* IERCOD=1 shows that the matrix is of lower rank than NORDRE */
8065 /* (the system is singular). */
8067 /* COMMONS USED : */
8068 /* ---------------- */
8070 /* REFERENCES CALLED : */
8071 /* ----------------------- */
8073 /* DESCRIPTION/NOTES/LIMITATIONS : */
8074 /* ----------------------------------- */
8075 /* ATTENTION : the indices of line and column are inverted */
8076 /* compared to usual indices. */
8078 /* a1*x + b1*y = c1 */
8079 /* a2*x + b2*y = c2 */
8080 /* should be represented by matrix ABMATR : */
8082 /* ABMATR(1,1) = a1 ABMATR(1,2) = a2 */
8083 /* ABMATR(2,1) = b1 ABMATR(2,2) = b2 */
8084 /* ABMATR(3,1) = c1 ABMATR(3,2) = c2 */
8086 /* To solve this system, it is necessary to set : */
8088 /* NORDRE = 2 (there are 2 equations with 2 unknown values), */
8089 /* NDIMEN = 1 (there is only one second member), */
8090 /* any NORMAX can be taken >= NORDRE. */
8092 /* To use this routine, it is recommended to use one of */
8093 /* interfaces : MMRSLWI or MMMRSLWD. */
8095 /* **********************************************************************
8098 /* Name of the routine */
8100 /* INTEGER IBB,MNFNDEB */
8103 /* IF (IBB.GE.2) CALL MGENMSG(NOMPR) */
8104 /* Parameter adjustments */
8105 xmatri_dim1 = *normax;
8106 xmatri_offset = xmatri_dim1 + 1;
8107 xmatri -= xmatri_offset;
8108 abmatr_dim1 = *nordre + *ndimen;
8109 abmatr_offset = abmatr_dim1 + 1;
8110 abmatr -= abmatr_offset;
8115 /* *********************************************************************
8117 /* Triangulation of matrix ABMATR. */
8118 /* *********************************************************************
8122 for (kk = 1; kk <= i__1; ++kk) {
8124 /* ---------- Find max pivot in column KK. ------------
8130 for (jj = kk; jj <= i__2; ++jj) {
8131 akj = (d__1 = abmatr[kk + jj * abmatr_dim1], advapp_abs(d__1));
8142 /* --------- Swapping of line KPIV with line KK. ------
8146 i__2 = *nordre + *ndimen;
8147 for (jj = kk; jj <= i__2; ++jj) {
8148 akj = abmatr[jj + kk * abmatr_dim1];
8149 abmatr[jj + kk * abmatr_dim1] = abmatr[jj + kpiv *
8151 abmatr[jj + kpiv * abmatr_dim1] = akj;
8156 /* ---------- Removal and triangularization. -----------
8159 pivot = -abmatr[kk + kk * abmatr_dim1];
8161 for (ii = kk + 1; ii <= i__2; ++ii) {
8162 akj = abmatr[kk + ii * abmatr_dim1] / pivot;
8163 i__3 = *nordre + *ndimen;
8164 for (jj = kk + 1; jj <= i__3; ++jj) {
8165 abmatr[jj + ii * abmatr_dim1] += akj * abmatr[jj + kk *
8176 /* *********************************************************************
8178 /* Solution of the system of triangular equations. */
8179 /* Matrix ABMATR(NORDRE+JJ,II), contains second members */
8180 /* of the system for 1<=j<=NDIMEN and 1<=i<=NORDRE. */
8181 /* *********************************************************************
8185 /* ---------------- Calculation of solutions by ascending. -----------------
8188 for (kk = *nordre; kk >= 1; --kk) {
8189 pivot = abmatr[kk + kk * abmatr_dim1];
8191 for (ii = 1; ii <= i__1; ++ii) {
8192 akj = abmatr[ii + *nordre + kk * abmatr_dim1];
8194 for (jj = kk + 1; jj <= i__2; ++jj) {
8195 akj -= abmatr[jj + kk * abmatr_dim1] * xmatri[jj + ii *
8199 xmatri[kk + ii * xmatri_dim1] = akj / pivot;
8206 /* ------If the absolute value of a pivot is smaller than --------
8207 /* ---------- EPSPIV: return the code of error. ------------
8217 AdvApp2Var_SysBase::maermsg_("MMRSLW ", iercod, 7L);
8219 /* IF (IBB.GE.2) CALL MGSOMSG(NOMPR) */
8223 //=======================================================================
8224 //function : AdvApp2Var_MathBase::mmmrslwd_
8226 //=======================================================================
8227 int AdvApp2Var_MathBase::mmmrslwd_(integer *normax,
8238 /* System generated locals */
8239 integer amat_dim1, amat_offset, bmat_dim1, bmat_offset, xmat_dim1,
8240 xmat_offset, aaux_dim1, aaux_offset, i__1, i__2;
8242 /* Local variables */
8243 static integer i__, j;
8246 /* IMPLICIT DOUBLE PRECISION (A-H,O-Z) */
8247 /* IMPLICIT INTEGER (I-N) */
8250 /* **********************************************************************
8255 /* Solution of a linear system by Gauss method where */
8256 /* the second member is a table of vectors. Method of partial pivot. */
8260 /* ALL, MATH_ACCES :: */
8261 /* SYSTEME&,EQUATION&, RESOLUTION,GAUSS ,&VECTEUR */
8263 /* INPUT ARGUMENTS : */
8264 /* ------------------ */
8265 /* NORMAX : Max. Dimension of AMAT. */
8266 /* NORDRE : Order of the matrix. */
8267 /* NDIM : Number of columns of BMAT and XMAT. */
8268 /* AMAT(NORMAX,NORDRE) : The processed matrix. */
8269 /* BMAT(NORMAX,NDIM) : The matrix of second member. */
8270 /* XMAT(NORMAX,NDIM) : The matrix of solutions. */
8271 /* EPSPIV : Min value of a pivot. */
8273 /* OUTPUT ARGUMENTS : */
8274 /* ------------------- */
8275 /* AAUX(NORDRE+NDIM,NORDRE) : Auxiliary matrix. */
8276 /* XMAT(NORMAX,NDIM) : Matrix of solutions. */
8277 /* IERCOD=0 shows that solutions in XMAT are valid. */
8278 /* IERCOD=1 shows that matrix AMAT is of lower rank than NORDRE. */
8280 /* COMMONS USED : */
8281 /* ---------------- */
8285 /* REFERENCES CALLED : */
8286 /* ---------------------- */
8288 /* MAERMSG MGENMSG MGSOMSG */
8289 /* MMRSLW I*4 MNFNDEB */
8291 /* DESCRIPTION/NOTES/LIMITATIONS : */
8292 /* ----------------------------------- */
8293 /* ATTENTION : lines and columns are located in usual order : */
8294 /* 1st index = index line */
8295 /* 2nd index = index column */
8296 /* Example, the system : */
8297 /* a1*x + b1*y = c1 */
8298 /* a2*x + b2*y = c2 */
8299 /* is represented by matrix AMAT : */
8301 /* AMAT(1,1) = a1 AMAT(2,1) = a2 */
8302 /* AMAT(1,2) = b1 AMAT(2,2) = b2 */
8304 /* The first index is the index of line, the second index */
8305 /* is the index of columns (Compare with MMRSLWI which is faster). */
8308 /* **********************************************************************
8311 /* Name of the routine */
8313 /* Parameter adjustments */
8314 amat_dim1 = *normax;
8315 amat_offset = amat_dim1 + 1;
8316 amat -= amat_offset;
8317 xmat_dim1 = *normax;
8318 xmat_offset = xmat_dim1 + 1;
8319 xmat -= xmat_offset;
8320 aaux_dim1 = *nordre + *ndim;
8321 aaux_offset = aaux_dim1 + 1;
8322 aaux -= aaux_offset;
8323 bmat_dim1 = *normax;
8324 bmat_offset = bmat_dim1 + 1;
8325 bmat -= bmat_offset;
8328 ibb = AdvApp2Var_SysBase::mnfndeb_();
8330 AdvApp2Var_SysBase::mgenmsg_("MMMRSLW", 7L);
8333 /* Initialization of the auxiliary matrix. */
8336 for (i__ = 1; i__ <= i__1; ++i__) {
8338 for (j = 1; j <= i__2; ++j) {
8339 aaux[j + i__ * aaux_dim1] = amat[i__ + j * amat_dim1];
8345 /* Second member. */
8348 for (i__ = 1; i__ <= i__1; ++i__) {
8350 for (j = 1; j <= i__2; ++j) {
8351 aaux[j + *nordre + i__ * aaux_dim1] = bmat[i__ + j * bmat_dim1];
8357 /* Solution of the system of equations. */
8359 mmrslw_(normax, nordre, ndim, epspiv, &aaux[aaux_offset], &xmat[
8360 xmat_offset], iercod);
8364 AdvApp2Var_SysBase::maermsg_("MMMRSLW", iercod, 7L);
8367 AdvApp2Var_SysBase::mgsomsg_("MMMRSLW", 7L);
8372 //=======================================================================
8373 //function : AdvApp2Var_MathBase::mmrtptt_
8375 //=======================================================================
8376 int AdvApp2Var_MathBase::mmrtptt_(integer *ndglgd,
8380 static integer ideb, nmod2, nsur2, ilong, ibb;
8383 /* **********************************************************************
8388 /* Extracts from Common LDGRTL the STRICTLY positive roots of the */
8389 /* Legendre polynom of degree NDGLGD, for 2 <= NDGLGD <= 61. */
8393 /* TOUS, AB_SPECIFI::COMMON&, EXTRACTION, &RACINE, &LEGENDRE. */
8395 /* INPUT ARGUMENTS : */
8396 /* ------------------ */
8397 /* NDGLGD : Mathematic degree of Legendre polynom. */
8398 /* This degree should be above or equal to 2 and */
8399 /* below or equal to 61. */
8401 /* OUTPUT ARGUMENTS : */
8402 /* ------------------- */
8403 /* RTLEGD : The table of strictly positive roots of */
8404 /* Legendre polynom of degree NDGLGD. */
8406 /* COMMONS USED : */
8407 /* ---------------- */
8409 /* REFERENCES CALLED : */
8410 /* ----------------------- */
8412 /* DESCRIPTION/NOTES/LIMITATIONS : */
8413 /* ----------------------------------- */
8414 /* ATTENTION: the condition on NDEGRE ( 2 <= NDEGRE <= 61) is not */
8415 /* tested. The caller should make the test. */
8418 /* **********************************************************************
8420 /* Nome of the routine */
8423 /* Common MLGDRTL: */
8424 /* This common includes POSITIVE roots of Legendre polynoms */
8425 /* AND the weight of Gauss quadrature formulas on all */
8426 /* POSITIVE roots of Legendre polynoms. */
8429 /* ***********************************************************************
8434 /* The common of Legendre roots. */
8440 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
8441 /* ----------------------------------- */
8444 /* ***********************************************************************
8450 /* ROOTAB : Table of all rotts of Legendre polynoms */
8451 /* between [0,1]. They are ranked for degrees increasing from 2 to 61. */
8452 /* HILTAB : Table of Legendre interpolators concerning ROOTAB. */
8453 /* The address is the same. */
8454 /* HI0TAB : Table of Legendre interpolators for root x=0 */
8455 /* the polynoms of UNEVEN degree. */
8456 /* RTLTB0 : Table of Li(uk) where uk are roots of a */
8457 /* Legendre polynom of EVEN degree. */
8458 /* RTLTB1 : Table of Li(uk) where uk are roots of a */
8459 /* Legendre polynom of UNEVEN degree. */
8462 /************************************************************************
8464 /* Parameter adjustments */
8468 ibb = AdvApp2Var_SysBase::mnfndeb_();
8470 AdvApp2Var_SysBase::mgenmsg_("MMRTPTT", 7L);
8476 nsur2 = *ndglgd / 2;
8477 nmod2 = *ndglgd % 2;
8480 ideb = nsur2 * (nsur2 - 1) / 2 + 1;
8481 AdvApp2Var_SysBase::mcrfill_((integer *)&ilong,
8482 (char *)&mlgdrtl_.rootab[ideb + nmod2 * 465 - 1],
8483 (char *)&rtlegd[1]);
8485 /* ----------------------------- The end --------------------------------
8490 AdvApp2Var_SysBase::mgsomsg_("MMRTPTT", 7L);
8495 //=======================================================================
8496 //function : AdvApp2Var_MathBase::mmsrre2_
8498 //=======================================================================
8499 int AdvApp2Var_MathBase::mmsrre2_(doublereal *tparam,
8507 /* System generated locals */
8510 /* Local variables */
8511 static integer ideb, ifin, imil, ibb;
8513 /* ***********************************************************************
8519 /* Find the interval corresponding to a valueb given in */
8520 /* increasing order of real numbers with double precision. */
8524 /* TOUS,MATH_ACCES::TABLEAU&,POINT&,CORRESPONDANCE,&RANG */
8526 /* INPUT ARGUMENTS : */
8527 /* ------------------ */
8529 /* TPARAM : Value to be tested. */
8530 /* NBRVAL : Size of TABLEV */
8531 /* TABLEV : Table of reals. */
8532 /* EPSIL : Epsilon of precision */
8534 /* OUTPUT ARGUMENTS : */
8535 /* ------------------- */
8537 /* NUMINT : Number of the interval (between 1 and NBRVAL-1). */
8538 /* ITYPEN : = 0 TPARAM is inside the interval NUMINT */
8539 /* = 1 : TPARAM corresponds to the lower limit of */
8540 /* the provided interval. */
8541 /* = 2 : TPARAM corresponds to the upper limit of */
8542 /* the provided interval. */
8544 /* IERCOD : Error code. */
8546 /* = 1 : TABLEV does not contain enough elements. */
8547 /* = 2 : TPARAM out of limits of TABLEV. */
8549 /* COMMONS USED : */
8550 /* ---------------- */
8552 /* REFERENCES CALLED : */
8553 /* ------------------- */
8555 /* DESCRIPTION/NOTES/LIMITATIONS : */
8556 /* --------------------------------- */
8557 /* There are NBRVAL values in TABLEV which stands for NBRVAL-1 intervals. */
8558 /* One searches the interval containing TPARAM by */
8559 /* dichotomy. Complexity of the algorithm : Log(n)/Log(2).(RBD). */
8561 /* ***********************************************************************
8565 /* Initialisations */
8567 /* Parameter adjustments */
8571 ibb = AdvApp2Var_SysBase::mnfndeb_();
8573 AdvApp2Var_SysBase::mgenmsg_("MMSRRE2", 7L);
8582 /* TABLEV should contain at least two values */
8589 /* TPARAM should be between extreme limits of TABLEV. */
8591 if (*tparam < tablev[1] || *tparam > tablev[*nbrval]) {
8596 /* ----------------------- SEARCH OF THE INTERVAL --------------------
8601 /* Test end of loop (found). */
8603 if (ideb + 1 == ifin) {
8608 /* Find by dichotomy on increasing values of TABLEV. */
8610 imil = (ideb + ifin) / 2;
8611 if (*tparam >= tablev[ideb] && *tparam <= tablev[imil]) {
8619 /* -------------- TEST IF TPARAM IS NOT A VALUE ---------
8620 /* ------------------------OF TABLEV UP TO EPSIL ----------------------
8624 if ((d__1 = *tparam - tablev[ideb], advapp_abs(d__1)) < *epsil) {
8628 if ((d__1 = *tparam - tablev[ifin], advapp_abs(d__1)) < *epsil) {
8633 /* --------------------------- THE END ----------------------------------
8638 AdvApp2Var_SysBase::maermsg_("MMSRRE2", iercod, 7L);
8641 AdvApp2Var_SysBase::mgsomsg_("MMSRRE2", 7L);
8646 //=======================================================================
8647 //function : mmtmave_
8649 //=======================================================================
8650 int mmtmave_(integer *nligne,
8660 /* System generated locals */
8663 /* Local variables */
8664 static logical ldbg;
8665 static integer imin, imax, i__, j, k;
8666 static doublereal somme;
8670 /* ***********************************************************************
8676 /* CREATES PRODUCT G V */
8677 /* WHERE THE MATRIX IS IN FORM OF PROFILE */
8681 /* RESERVE, PRODUCT, MATRIX, PROFILE, VECTOR */
8683 /* INPUT ARGUMENTS : */
8684 /* -------------------- */
8685 /* NLIGNE : NUMBER OF LINE OF THE MATRIX */
8686 /* NCOLON : NOMBER OF COLUMN OF THE MATRIX */
8687 /* GPOSIT: TABLE OF POSITIONING OF TERMS OF STORAGE */
8688 /* GPOSIT(1,I) CONTAINS THE NUMBER of TERMS-1 ON LINE
8689 I IN THE PROFILE OF THE MATRIX */
8690 /* GPOSIT(2,I) CONTAINS THE INDEX OF STORAGE OF THE DIAGONAL TERM
8692 /* GPOSIT(3,I) CONTAINS THE INDEX COLUMN OF THE FIRST TERM OF
8693 /* PROFILE OF LINE I */
8694 /* GNSTOC : NOMBER OF TERM IN THE PROFILE OF GMATRI */
8695 /* GMATRI : MATRIX OF CONSTRAINTS IN FORM OF PROFILE */
8696 /* VECIN : INPUT VECTOR */
8698 /* OUTPUT ARGUMENTS : */
8699 /* --------------------- */
8700 /* VECOUT : VECTOR PRODUCT */
8701 /* IERCOD : ERROR CODE */
8704 /* COMMONS USED : */
8705 /* ------------------ */
8708 /* REFERENCES CALLED : */
8709 /* --------------------- */
8712 /* DESCRIPTION/NOTES/LIMITATIONS : */
8713 /* ----------------------------------- */
8715 /* ***********************************************************************
8718 /* ***********************************************************************
8723 /* ***********************************************************************
8725 /* INITIALISATIONS */
8726 /* ***********************************************************************
8729 /* Parameter adjustments */
8736 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
8738 AdvApp2Var_SysBase::mgenmsg_("MMTMAVE", 7L);
8742 /* ***********************************************************************
8745 /* ***********************************************************************
8751 for (i__ = 1; i__ <= i__1; ++i__) {
8754 for (j = 1; j <= i__2; ++j) {
8755 imin = gposit[j * 3 + 3];
8756 imax = gposit[j * 3 + 1] + gposit[j * 3 + 3] - 1;
8757 aux = gposit[j * 3 + 2] - gposit[j * 3 + 1] - imin + 1;
8758 if (imin <= i__ && i__ <= imax) {
8760 somme += gmatri[k] * vecin[j];
8763 vecout[i__] = somme;
8772 /* ***********************************************************************
8774 /* ERROR PROCESSING */
8775 /* ***********************************************************************
8779 /* ***********************************************************************
8781 /* RETURN CALLING PROGRAM */
8782 /* ***********************************************************************
8787 /* ___ DESALLOCATION, ... */
8789 AdvApp2Var_SysBase::maermsg_("MMTMAVE", iercod, 7L);
8791 AdvApp2Var_SysBase::mgsomsg_("MMTMAVE", 7L);
8796 //=======================================================================
8797 //function : mmtrpj0_
8799 //=======================================================================
8800 int mmtrpj0_(integer *ncofmx,
8810 /* System generated locals */
8811 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
8814 /* Local variables */
8815 static integer ncut, i__;
8816 static doublereal bidon, error;
8820 /* ***********************************************************************
8825 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
8826 /* Legendre with a given precision. */
8830 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
8832 /* INPUT ARGUMENTS : */
8833 /* ------------------ */
8834 /* NCOFMX : Max Nb of coeff. of the curve (dimensioning). */
8835 /* NDIMEN : Dimension of the space. */
8836 /* NCOEFF : Degree +1 of the polynom. */
8837 /* EPSI3D : Precision required for the approximation. */
8838 /* CRVLGD : The curve the degree which of it is required to lower. */
8840 /* OUTPUT ARGUMENTS : */
8841 /* ------------------- */
8842 /* EPSTRC : Precision of the approximation. */
8843 /* NCFNEW : Degree +1 of the resulting polynom. */
8845 /* COMMONS USED : */
8846 /* ---------------- */
8848 /* REFERENCES CALLED : */
8849 /* ----------------------- */
8851 /* DESCRIPTION/NOTES/LIMITATIONS : */
8852 /* ----------------------------------- */
8854 /* ***********************************************************************
8858 /* ------- Minimum degree that can be attained : Stop at 1 (RBD) ---------
8861 /* Parameter adjustments */
8863 crvlgd_dim1 = *ncofmx;
8864 crvlgd_offset = crvlgd_dim1 + 1;
8865 crvlgd -= crvlgd_offset;
8869 /* ------------------- Init for error calculation -----------------------
8872 for (i__ = 1; i__ <= i__1; ++i__) {
8879 /* Cutting of coefficients. */
8882 /* ------ Loop on the series of Legendre :NCOEFF --> 2 (RBD) -----------
8885 for (i__ = *ncoeff; i__ >= i__1; --i__) {
8886 /* Factor of renormalization. */
8887 bidon = ((i__ - 1) * 2. + 1.) / 2.;
8888 bidon = sqrt(bidon);
8890 for (nd = 1; nd <= i__2; ++nd) {
8891 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
8895 /* Cutting is stopped if the norm becomes too great. */
8896 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
8897 if (error > *epsi3d) {
8902 /* --- Max error cumulee when the I-th coeff is removed. */
8909 /* --------------------------------- End --------------------------------
8916 //=======================================================================
8917 //function : mmtrpj2_
8919 //=======================================================================
8920 int mmtrpj2_(integer *ncofmx,
8930 /* Initialized data */
8932 static doublereal xmaxj[57] = { .9682458365518542212948163499456,
8933 .986013297183269340427888048593603,
8934 1.07810420343739860362585159028115,
8935 1.17325804490920057010925920756025,
8936 1.26476561266905634732910520370741,
8937 1.35169950227289626684434056681946,
8938 1.43424378958284137759129885012494,
8939 1.51281316274895465689402798226634,
8940 1.5878364329591908800533936587012,
8941 1.65970112228228167018443636171226,
8942 1.72874345388622461848433443013543,
8943 1.7952515611463877544077632304216,
8944 1.85947199025328260370244491818047,
8945 1.92161634324190018916351663207101,
8946 1.98186713586472025397859895825157,
8947 2.04038269834980146276967984252188,
8948 2.09730119173852573441223706382076,
8949 2.15274387655763462685970799663412,
8950 2.20681777186342079455059961912859,
8951 2.25961782459354604684402726624239,
8952 2.31122868752403808176824020121524,
8953 2.36172618435386566570998793688131,
8954 2.41117852396114589446497298177554,
8955 2.45964731268663657873849811095449,
8956 2.50718840313973523778244737914028,
8957 2.55385260994795361951813645784034,
8958 2.59968631659221867834697883938297,
8959 2.64473199258285846332860663371298,
8960 2.68902863641518586789566216064557,
8961 2.73261215675199397407027673053895,
8962 2.77551570192374483822124304745691,
8963 2.8177699459714315371037628127545,
8964 2.85940333797200948896046563785957,
8965 2.90044232019793636101516293333324,
8966 2.94091151970640874812265419871976,
8967 2.98083391718088702956696303389061,
8968 3.02023099621926980436221568258656,
8969 3.05912287574998661724731962377847,
8970 3.09752842783622025614245706196447,
8971 3.13546538278134559341444834866301,
8972 3.17295042316122606504398054547289,
8973 3.2099992681699613513775259670214,
8974 3.24662674946606137764916854570219,
8975 3.28284687953866689817670991319787,
8976 3.31867291347259485044591136879087,
8977 3.35411740487202127264475726990106,
8978 3.38919225660177218727305224515862,
8979 3.42390876691942143189170489271753,
8980 3.45827767149820230182596660024454,
8981 3.49230918177808483937957161007792,
8982 3.5260130200285724149540352829756,
8983 3.55939845146044235497103883695448,
8984 3.59247431368364585025958062194665,
8985 3.62524904377393592090180712976368,
8986 3.65773070318071087226169680450936,
8987 3.68992700068237648299565823810245,
8988 3.72184531357268220291630708234186 };
8990 /* System generated locals */
8991 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
8994 /* Local variables */
8995 static integer ncut, i__;
8996 static doublereal bidon, error;
8997 static integer ia, nd;
8998 static doublereal bid, eps1;
9001 /* ***********************************************************************
9006 /* Lower the degree of a curve defined on (-1,1) in the direction of */
9007 /* Legendre with a given precision. */
9011 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
9013 /* INPUT ARGUMENTS : */
9014 /* ------------------ */
9015 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9016 /* NDIMEN : Dimension of the space. */
9017 /* NCOEFF : Degree +1 of the polynom. */
9018 /* EPSI3D : Precision required for the approximation. */
9019 /* CRVLGD : The curve the degree which of will be lowered. */
9021 /* OUTPUT ARGUMENTS : */
9022 /* ------------------- */
9023 /* YCVMAX : Auxiliary table (error max on each dimension).
9025 /* EPSTRC : Precision of the approximation. */
9026 /* NCFNEW : Degree +1 of the resulting polynom. */
9028 /* COMMONS USED : */
9029 /* ---------------- */
9031 /* REFERENCES CALLED : */
9032 /* ----------------------- */
9034 /* DESCRIPTION/NOTES/LIMITATIONS : */
9035 /* ----------------------------------- */
9037 /* ***********************************************************************
9041 /* Parameter adjustments */
9043 crvlgd_dim1 = *ncofmx;
9044 crvlgd_offset = crvlgd_dim1 + 1;
9045 crvlgd -= crvlgd_offset;
9051 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9055 /* Init for calculation of error. */
9057 for (i__ = 1; i__ <= i__1; ++i__) {
9064 /* Cutting of coefficients. */
9067 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9070 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9071 /* Factor of renormalization. */
9072 bidon = xmaxj[i__ - ncut];
9074 for (nd = 1; nd <= i__2; ++nd) {
9075 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
9079 /* One stops to cut if the norm becomes too great. */
9080 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9081 if (error > *epsi3d) {
9086 /* --- Max error cumulated when the I-th coeff is removed. */
9093 /* ------- Cutting of zero coeffs of interpolation (RBD) -------
9097 if (*ncfnew == ia) {
9098 AdvApp2Var_MathBase::mmeps1_(&eps1);
9099 for (i__ = ia; i__ >= 2; --i__) {
9102 for (nd = 1; nd <= i__1; ++nd) {
9103 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
9112 /* --- If all coeffs can be removed, this is a point. */
9116 /* --------------------------------- End --------------------------------
9123 //=======================================================================
9124 //function : mmtrpj4_
9126 //=======================================================================
9127 int mmtrpj4_(integer *ncofmx,
9136 /* Initialized data */
9138 static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
9139 1.05299572648705464724876659688996,
9140 1.0949715351434178709281698645813,
9141 1.15078388379719068145021100764647,
9142 1.2094863084718701596278219811869,
9143 1.26806623151369531323304177532868,
9144 1.32549784426476978866302826176202,
9145 1.38142537365039019558329304432581,
9146 1.43575531950773585146867625840552,
9147 1.48850442653629641402403231015299,
9148 1.53973611681876234549146350844736,
9149 1.58953193485272191557448229046492,
9150 1.63797820416306624705258190017418,
9151 1.68515974143594899185621942934906,
9152 1.73115699602477936547107755854868,
9153 1.77604489805513552087086912113251,
9154 1.81989256661534438347398400420601,
9155 1.86276344480103110090865609776681,
9156 1.90471563564740808542244678597105,
9157 1.94580231994751044968731427898046,
9158 1.98607219357764450634552790950067,
9159 2.02556989246317857340333585562678,
9160 2.06433638992049685189059517340452,
9161 2.10240936014742726236706004607473,
9162 2.13982350649113222745523925190532,
9163 2.17661085564771614285379929798896,
9164 2.21280102016879766322589373557048,
9165 2.2484214321456956597803794333791,
9166 2.28349755104077956674135810027654,
9167 2.31805304852593774867640120860446,
9168 2.35210997297725685169643559615022,
9169 2.38568889602346315560143377261814,
9170 2.41880904328694215730192284109322,
9171 2.45148841120796359750021227795539,
9172 2.48374387161372199992570528025315,
9173 2.5155912654873773953959098501893,
9174 2.54704548720896557684101746505398,
9175 2.57812056037881628390134077704127,
9176 2.60882970619319538196517982945269,
9177 2.63918540521920497868347679257107,
9178 2.66919945330942891495458446613851,
9179 2.69888301230439621709803756505788,
9180 2.72824665609081486737132853370048,
9181 2.75730041251405791603760003778285,
9182 2.78605380158311346185098508516203,
9183 2.81451587035387403267676338931454,
9184 2.84269522483114290814009184272637,
9185 2.87060005919012917988363332454033,
9186 2.89823818258367657739520912946934,
9187 2.92561704377132528239806135133273,
9188 2.95274375377994262301217318010209,
9189 2.97962510678256471794289060402033,
9190 3.00626759936182712291041810228171,
9191 3.03267744830655121818899164295959,
9192 3.05886060707437081434964933864149 };
9194 /* System generated locals */
9195 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
9198 /* Local variables */
9199 static integer ncut, i__;
9200 static doublereal bidon, error;
9201 static integer ia, nd;
9202 static doublereal bid, eps1;
9206 /* ***********************************************************************
9211 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
9212 /* Legendre with a given precision. */
9216 /* LEGENDRE, POLYGON, TRONCATION, CURVE, SMOOTHING. */
9218 /* INPUT ARGUMENTS : */
9219 /* ------------------ */
9220 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9221 /* NDIMEN : Dimension of the space. */
9222 /* NCOEFF : Degree +1 of the polynom. */
9223 /* EPSI3D : Precision required for the approximation. */
9224 /* CRVLGD : The curve which wishes to lower the degree. */
9226 /* OUTPUT ARGUMENTS : */
9227 /* ------------------- */
9228 /* YCVMAX : Auxiliary table (max error on each dimension).
9230 /* EPSTRC : Precision of the approximation. */
9231 /* NCFNEW : Degree +1 of the resulting polynom. */
9233 /* COMMONS USED : */
9234 /* ---------------- */
9236 /* REFERENCES CALLED : */
9237 /* ----------------------- */
9239 /* DESCRIPTION/NOTES/LIMITATIONS : */
9240 /* ----------------------------------- */
9242 /* ***********************************************************************
9246 /* Parameter adjustments */
9248 crvlgd_dim1 = *ncofmx;
9249 crvlgd_offset = crvlgd_dim1 + 1;
9250 crvlgd -= crvlgd_offset;
9256 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9260 /* Init for error calculation. */
9262 for (i__ = 1; i__ <= i__1; ++i__) {
9269 /* Cutting of coefficients. */
9272 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9275 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9276 /* Factor of renormalization. */
9277 bidon = xmaxj[i__ - ncut];
9279 for (nd = 1; nd <= i__2; ++nd) {
9280 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
9284 /* Stop cutting if the norm becomes too great. */
9285 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9286 if (error > *epsi3d) {
9291 /* -- Error max cumulated when the I-eme coeff is removed. */
9298 /* ------- Cutting of zero coeffs of the pole of interpolation (RBD) -------
9302 if (*ncfnew == ia) {
9303 AdvApp2Var_MathBase::mmeps1_(&eps1);
9304 for (i__ = ia; i__ >= 2; --i__) {
9307 for (nd = 1; nd <= i__1; ++nd) {
9308 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
9317 /* --- If all coeffs can be removed, this is a point. */
9321 /* --------------------------------- End --------------------------------
9328 //=======================================================================
9329 //function : mmtrpj6_
9331 //=======================================================================
9332 int mmtrpj6_(integer *ncofmx,
9342 /* Initialized data */
9344 static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
9345 1.11626917091567929907256116528817,
9346 1.1327140810290884106278510474203,
9347 1.1679452722668028753522098022171,
9348 1.20910611986279066645602153641334,
9349 1.25228283758701572089625983127043,
9350 1.29591971597287895911380446311508,
9351 1.3393138157481884258308028584917,
9352 1.3821288728999671920677617491385,
9353 1.42420414683357356104823573391816,
9354 1.46546895108549501306970087318319,
9355 1.50590085198398789708599726315869,
9356 1.54550385142820987194251585145013,
9357 1.58429644271680300005206185490937,
9358 1.62230484071440103826322971668038,
9359 1.65955905239130512405565733793667,
9360 1.69609056468292429853775667485212,
9361 1.73193098017228915881592458573809,
9362 1.7671112206990325429863426635397,
9363 1.80166107681586964987277458875667,
9364 1.83560897003644959204940535551721,
9365 1.86898184653271388435058371983316,
9366 1.90180515174518670797686768515502,
9367 1.93410285411785808749237200054739,
9368 1.96589749778987993293150856865539,
9369 1.99721027139062501070081653790635,
9370 2.02806108474738744005306947877164,
9371 2.05846864831762572089033752595401,
9372 2.08845055210580131460156962214748,
9373 2.11802334209486194329576724042253,
9374 2.14720259305166593214642386780469,
9375 2.17600297710595096918495785742803,
9376 2.20443832785205516555772788192013,
9377 2.2325216999457379530416998244706,
9378 2.2602654243075083168599953074345,
9379 2.28768115912702794202525264301585,
9380 2.3147799369092684021274946755348,
9381 2.34157220782483457076721300512406,
9382 2.36806787963276257263034969490066,
9383 2.39427635443992520016789041085844,
9384 2.42020656255081863955040620243062,
9385 2.44586699364757383088888037359254,
9386 2.47126572552427660024678584642791,
9387 2.49641045058324178349347438430311,
9388 2.52130850028451113942299097584818,
9389 2.54596686772399937214920135190177,
9390 2.5703922285006754089328998222275,
9391 2.59459096001908861492582631591134,
9392 2.61856915936049852435394597597773,
9393 2.64233265984385295286445444361827,
9394 2.66588704638685848486056711408168,
9395 2.68923766976735295746679957665724,
9396 2.71238965987606292679677228666411 };
9398 /* System generated locals */
9399 integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
9402 /* Local variables */
9403 static integer ncut, i__;
9404 static doublereal bidon, error;
9405 static integer ia, nd;
9406 static doublereal bid, eps1;
9410 /* ***********************************************************************
9415 /* Lowers the degree of a curve defined on (-1,1) in the direction of */
9416 /* Legendre to a given precision. */
9420 /* LEGENDRE,POLYGON,TRUNCATION,CURVE,SMOOTHING. */
9422 /* INPUT ARGUMENTS : */
9423 /* ------------------ */
9424 /* NCOFMX : Max nb of coeff. of the curve (dimensioning). */
9425 /* NDIMEN : Dimension of the space. */
9426 /* NCOEFF : Degree +1 of the polynom. */
9427 /* EPSI3D : Precision required for the approximation. */
9428 /* CRVLGD : The curve the degree which of will be lowered. */
9430 /* OUTPUT ARGUMENTS : */
9431 /* ------------------- */
9432 /* YCVMAX : Auxiliary table (max error on each dimension).
9433 /* EPSTRC : Precision of the approximation. */
9434 /* NCFNEW : Degree +1 of the resulting polynom. */
9436 /* COMMONS USED : */
9437 /* ---------------- */
9439 /* REFERENCES CALLED : */
9440 /* ----------------------- */
9442 /* DESCRIPTION/NOTES/LIMITATIONS : */
9443 /* ----------------------------------- */
9445 /* ***********************************************************************
9449 /* Parameter adjustments */
9451 crvlgd_dim1 = *ncofmx;
9452 crvlgd_offset = crvlgd_dim1 + 1;
9453 crvlgd -= crvlgd_offset;
9459 /* Minimum degree that can be reached : Stop at IA (RBD). -------------
9463 /* Init for error calculation. */
9465 for (i__ = 1; i__ <= i__1; ++i__) {
9472 /* Cutting of coefficients. */
9475 /* ------ Loop on the series of Jacobi :NCOEFF --> IA+1 (RBD) ----------
9478 for (i__ = *ncoeff; i__ >= i__1; --i__) {
9479 /* Factor of renormalization. */
9480 bidon = xmaxj[i__ - ncut];
9482 for (nd = 1; nd <= i__2; ++nd) {
9483 ycvmax[nd] += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1)) *
9487 /* Stop cutting if the norm becomes too great. */
9488 error = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
9489 if (error > *epsi3d) {
9494 /* --- Max error cumulated when the I-th coeff is removed. */
9501 /* ------- Cutting of zero coeff. of the pole of interpolation (RBD) -------
9505 if (*ncfnew == ia) {
9506 AdvApp2Var_MathBase::mmeps1_(&eps1);
9507 for (i__ = ia; i__ >= 2; --i__) {
9510 for (nd = 1; nd <= i__1; ++nd) {
9511 bid += (d__1 = crvlgd[i__ + nd * crvlgd_dim1], advapp_abs(d__1));
9520 /* --- If all coeffs can be removed, this is a point. */
9524 /* --------------------------------- End --------------------------------
9531 //=======================================================================
9532 //function : AdvApp2Var_MathBase::mmtrpjj_
9534 //=======================================================================
9535 int AdvApp2Var_MathBase::mmtrpjj_(integer *ncofmx,
9545 /* System generated locals */
9546 integer crvlgd_dim1, crvlgd_offset;
9548 /* Local variables */
9552 /* ***********************************************************************
9557 /* Lower the degree of a curve defined on (-1,1) in the direction of */
9558 /* Legendre with a given precision. */
9562 /* LEGENDRE, POLYGON, TRUNCATION, CURVE, SMOOTHING. */
9564 /* INPUT ARGUMENTS : */
9565 /* ------------------ */
9566 /* NCOFMX : Max Nb coeff. of the curve (dimensioning). */
9567 /* NDIMEN : Dimension of the space. */
9568 /* NCOEFF : Degree +1 of the polynom. */
9569 /* EPSI3D : Precision required for the approximation. */
9570 /* IORDRE : Order of continuity at the extremities. */
9571 /* CRVLGD : The curve the degree which of should be lowered. */
9573 /* OUTPUT ARGUMENTS : */
9574 /* ------------------- */
9575 /* ERRMAX : Precision of the approximation. */
9576 /* NCFNEW : Degree +1 of the resulting polynom. */
9578 /* COMMONS USED : */
9579 /* ---------------- */
9581 /* REFERENCES CALLED : */
9582 /* ----------------------- */
9584 /* DESCRIPTION/NOTES/LIMITATIONS : */
9585 /* ----------------------------------- */
9587 /* ***********************************************************************
9591 /* Parameter adjustments */
9593 crvlgd_dim1 = *ncofmx;
9594 crvlgd_offset = crvlgd_dim1 + 1;
9595 crvlgd -= crvlgd_offset;
9598 ia = (*iordre + 1) << 1;
9601 mmtrpj0_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9602 ycvmax[1], errmax, ncfnew);
9603 } else if (ia == 2) {
9604 mmtrpj2_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9605 ycvmax[1], errmax, ncfnew);
9606 } else if (ia == 4) {
9607 mmtrpj4_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9608 ycvmax[1], errmax, ncfnew);
9610 mmtrpj6_(ncofmx, ndimen, ncoeff, epsi3d, &crvlgd[crvlgd_offset], &
9611 ycvmax[1], errmax, ncfnew);
9614 /* ------------------------ End -----------------------------------------
9620 //=======================================================================
9621 //function : AdvApp2Var_MathBase::mmunivt_
9623 //=======================================================================
9624 int AdvApp2Var_MathBase::mmunivt_(integer *ndimen,
9631 static doublereal c_b2 = 10.;
9633 /* System generated locals */
9637 /* Local variables */
9638 static integer nchif, iunit, izero;
9639 static doublereal vnorm;
9641 static doublereal bid;
9642 static doublereal eps0;
9647 /* ***********************************************************************
9652 /* CALCULATE THE NORMAL VECTOR BASING ON ANY VECTOR */
9653 /* WITH PRECISION GIVEN BY THE USER. */
9657 /* ALL, MATH_ACCES :: */
9658 /* VECTEUR&, NORMALISATION, &VECTEUR */
9660 /* INPUT ARGUMENTS : */
9661 /* ------------------ */
9662 /* NDIMEN : DIMENSION OF THE SPACE */
9663 /* VECTOR : VECTOR TO BE NORMED */
9664 /* EPSILN : EPSILON BELOW WHICH IT IS CONSIDERED THAT THE */
9665 /* NORM OF THE VECTOR IS NULL. IF EPSILN<=0, A DEFAULT VALUE */
9666 /* IS IMPOSED (10.D-17 ON VAX). */
9668 /* OUTPUT ARGUMENTS : */
9669 /* ------------------- */
9670 /* VECNRM : NORMED VECTOR */
9671 /* IERCOD 101 : THE VECTOR IS NULL UP TO EPSILN. */
9674 /* COMMONS USED : */
9675 /* ---------------- */
9677 /* REFERENCES CALLED : */
9678 /* ----------------------- */
9680 /* DESCRIPTION/NOTES/LIMITATIONS : */
9681 /* ----------------------------------- */
9682 /* VECTOR and VECNRM can be identic. */
9684 /* The norm of vector is calculated and each component is divided by
9685 /* this norm. After this it is checked if all componentes of the */
9686 /* vector except for one cost 0 with machine precision. In */
9687 /* this case the quasi-null components are set to 0.D0. */
9689 /* ***********************************************************************
9693 /* Parameter adjustments */
9700 /* -------- Precision by default : zero machine 10.D-17 on Vax ------
9703 AdvApp2Var_SysBase::maovsr8_(&nchif);
9704 if (*epsiln <= 0.) {
9706 eps0 = AdvApp2Var_MathBase::pow__di(&c_b2, &i__1);
9711 /* ------------------------- Calculation of the norm --------------------
9714 vnorm = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vector[1]);
9715 if (vnorm <= eps0) {
9716 AdvApp2Var_SysBase::mvriraz_((integer *)ndimen, (char *)&vecnrm[1]);
9721 /* ---------------------- Calculation of the vector norm ---------------
9725 i__1 = (-nchif - 1) / 2;
9726 eps0 = AdvApp2Var_MathBase::pow__di(&c_b2, &i__1);
9728 for (ii = 1; ii <= i__1; ++ii) {
9729 vecnrm[ii] = vector[ii] / vnorm;
9730 if ((d__1 = vecnrm[ii], advapp_abs(d__1)) <= eps0) {
9738 /* ------ Case when all coordinates except for one are almost null ----
9740 /* ------------- then one of coordinates costs 1.D0 or -1.D0 --------
9743 if (izero == *ndimen - 1) {
9744 bid = vecnrm[iunit];
9746 for (ii = 1; ii <= i__1; ++ii) {
9753 vecnrm[iunit] = -1.;
9757 /* -------------------------------- The end -----------------------------
9764 //=======================================================================
9765 //function : AdvApp2Var_MathBase::mmveps3_
9767 //=======================================================================
9768 int AdvApp2Var_MathBase::mmveps3_(doublereal *eps03)
9770 /* Initialized data */
9772 static char nomprg[8+1] = "MMEPS1 ";
9778 /************************************************************************
9783 /* Extraction of EPS1 from COMMON MPRCSN. */
9787 /* MPRCSN,PRECISON,EPS3. */
9789 /* INPUT ARGUMENTS : */
9790 /* ------------------ */
9793 /* OUTPUT ARGUMENTS : */
9794 /* ------------------- */
9795 /* EPS3 : space zero of the denominator (10**-9) */
9796 /* EPS3 should value 10**-15 */
9798 /* COMMONS USED : */
9799 /* ---------------- */
9801 /* REFERENCES CALLED : */
9802 /* ----------------------- */
9804 /* DESCRIPTION/NOTES/LIMITATIONS : */
9805 /* ----------------------------------- */
9808 /* ***********************************************************************
9813 /* ***********************************************************************
9818 /* GIVES TOLERANCES OF NULLITY IN STRIM */
9819 /* AND LIMITS OF ITERATIVE PROCESSES */
9821 /* GENERAL CONTEXT, MODIFIABLE BY THE UTILISER */
9825 /* PARAMETER , TOLERANCE */
9827 /* DESCRIPTION/NOTES/LIMITATIONS : */
9828 /* ----------------------------------- */
9829 /* INITIALISATION : PROFILE , **VIA MPRFTX** AT INPUT IN STRIM*/
9830 /* LOADING OF DEFAULT VALUES OF THE PROFILE IN MPRFTX AT INPUT*/
9831 /* IN STRIM. THEY ARE PRESERVED IN THE LOCAL VARIABLES OF MPRFTX */
9833 /* RESET DEFAULT VALUES : MDFINT */
9834 /* MODIFICATION INTERACTIVE BY THE USER : MDBINT */
9836 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
9837 /* MEPSPB ... EPS3,EPS4 */
9838 /* MEPSLN ... EPS2, NITERM , NITERR */
9839 /* MEPSNR ... EPS2 , NITERM */
9840 /* MITERR ... NITERR */
9843 /* ***********************************************************************
9846 /* NITERM : MAX NB OF ITERATIONS */
9847 /* NITERR : NB OF RAPID ITERATIONS */
9848 /* EPS1 : TOLERANCE OF 3D NULL DISTANCE */
9849 /* EPS2 : TOLERANCE OF ZERO PARAMETRIC DISTANCE */
9850 /* EPS3 : TOLERANCE TO AVOID DIVISION BY 0.. */
9851 /* EPS4 : TOLERANCE ANGULAR */
9855 /* ***********************************************************************
9858 ibb = AdvApp2Var_SysBase::mnfndeb_();
9860 AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
9863 *eps03 = mmprcsn_.eps3;
9868 //=======================================================================
9869 //function : AdvApp2Var_MathBase::mmvncol_
9871 //=======================================================================
9872 int AdvApp2Var_MathBase::mmvncol_(integer *ndimen,
9878 /* System generated locals */
9881 /* Local variables */
9882 static logical ldbg;
9884 static doublereal vaux1[3], vaux2[3];
9885 static logical colin;
9886 static doublereal valaux;
9890 /* ***********************************************************************
9895 /* CALCULATE A VECTOR NON-COLINEAR TO A GIVEN NON-NULL VECTOR */
9899 /* PUBLIC, VECTOR, FREE */
9901 /* INPUT ARGUMENTS : */
9902 /* -------------------- */
9903 /* ndimen : dimension of the space */
9904 /* vecin : input vector */
9906 /* OUTPUT ARGUMENTS : */
9907 /* --------------------- */
9909 /* vecout : vector non colinear to vecin */
9911 /* COMMONS USED : */
9912 /* ------------------ */
9915 /* REFERENCES CALLED : */
9916 /* --------------------- */
9919 /* DESCRIPTION/NOTES/LIMITATIONS : */
9920 /* ----------------------------------- */
9922 /* ***********************************************************************
9925 /* ***********************************************************************
9930 /* ***********************************************************************
9932 /* INITIALISATIONS */
9933 /* ***********************************************************************
9936 /* Parameter adjustments */
9941 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
9943 AdvApp2Var_SysBase::mgenmsg_("MMVNCOL", 7L);
9947 /* ***********************************************************************
9950 /* ***********************************************************************
9953 if (*ndimen <= 1 || *ndimen > 3) {
9959 while(d__ <= *ndimen) {
9960 if (vecin[d__] == 0.) {
9965 if (aux == *ndimen) {
9970 for (d__ = 1; d__ <= 3; ++d__) {
9971 vaux1[d__ - 1] = 0.;
9974 for (d__ = 1; d__ <= i__1; ++d__) {
9975 vaux1[d__ - 1] = vecin[d__];
9976 vaux2[d__ - 1] = vecin[d__];
9985 vaux2[d__ - 1] += 1;
9986 valaux = vaux1[1] * vaux2[2] - vaux1[2] * vaux2[1];
9988 valaux = vaux1[2] * vaux2[0] - vaux1[0] * vaux2[2];
9990 valaux = vaux1[0] * vaux2[1] - vaux1[1] * vaux2[0];
10005 for (d__ = 1; d__ <= i__1; ++d__) {
10006 vecout[d__] = vaux2[d__ - 1];
10011 /* ***********************************************************************
10013 /* ERROR PROCESSING */
10014 /* ***********************************************************************
10023 /* ***********************************************************************
10025 /* RETURN CALLING PROGRAM */
10026 /* ***********************************************************************
10032 AdvApp2Var_SysBase::maermsg_("MMVNCOL", iercod, 7L);
10034 AdvApp2Var_SysBase::mgsomsg_("MMVNCOL", 7L);
10039 //=======================================================================
10040 //function : AdvApp2Var_MathBase::mmwprcs_
10042 //=======================================================================
10043 void AdvApp2Var_MathBase::mmwprcs_(doublereal *epsil1,
10044 doublereal *epsil2,
10045 doublereal *epsil3,
10046 doublereal *epsil4,
10053 /* ***********************************************************************
10058 /* ACCESS IN WRITING FOR COMMON MPRCSN */
10064 /* INPUT ARGUMENTS : */
10065 /* -------------------- */
10066 /* EPSIL1 : TOLERANCE OF 3D NULL DISTANCE */
10067 /* EPSIL2 : TOLERANCE OF PARAMETRIC NULL DISTANCE */
10068 /* EPSIL3 : TOLERANCE TO AVOID DIVISION BY 0.. */
10069 /* EPSIL4 : ANGULAR TOLERANCE */
10070 /* NITER1 : MAX NB OF ITERATIONS */
10071 /* NITER2 : NB OF RAPID ITERATIONS */
10073 /* OUTPUT ARGUMENTS : */
10074 /* --------------------- */
10077 /* COMMONS USED : */
10078 /* ------------------ */
10081 /* REFERENCES CALLED : */
10082 /* --------------------- */
10085 /* DESCRIPTION/NOTES/LIMITATIONS : */
10086 /* ----------------------------------- */
10089 /* ***********************************************************************
10092 /* ***********************************************************************
10096 /* ***********************************************************************
10098 /* INITIALIZATIONS */
10099 /* ***********************************************************************
10102 /* ***********************************************************************
10105 /* ***********************************************************************
10108 /* ***********************************************************************
10113 /* GIVES TOLERANCES OF NULLITY IN STRIM */
10114 /* AND LIMITS OF ITERATIVE PROCESSES */
10116 /* GENERAL CONTEXT, MODIFIABLE BY THE UTILISER */
10120 /* PARAMETER , TOLERANCE */
10122 /* DESCRIPTION/NOTES/LIMITATIONS : */
10123 /* ----------------------------------- */
10124 /* INITIALISATION : PROFILE , **VIA MPRFTX** AT INPUT IN STRIM*/
10125 /* LOADING OF DEFAULT VALUES OF THE PROFILE IN MPRFTX AT INPUT*/
10126 /* IN STRIM. THEY ARE PRESERVED IN THE LOCAL VARIABLES OF MPRFTX */
10128 /* RESET DEFAULT VALUES : MDFINT */
10129 /* MODIFICATION INTERACTIVE BY THE USER : MDBINT */
10131 /* ACCESS FUNCTION : MMEPS1 ... EPS1 */
10132 /* MEPSPB ... EPS3,EPS4 */
10133 /* MEPSLN ... EPS2, NITERM , NITERR */
10134 /* MEPSNR ... EPS2 , NITERM */
10135 /* MITERR ... NITERR */
10138 /* ***********************************************************************
10141 /* NITERM : MAX NB OF ITERATIONS */
10142 /* NITERR : NB OF RAPID ITERATIONS */
10143 /* EPS1 : TOLERANCE OF 3D NULL DISTANCE */
10144 /* EPS2 : TOLERANCE OF ZERO PARAMETRIC DISTANCE */
10145 /* EPS3 : TOLERANCE TO AVOID DIVISION BY 0.. */
10146 /* EPS4 : TOLERANCE ANGULAR */
10149 /* ***********************************************************************
10151 mmprcsn_.eps1 = *epsil1;
10152 mmprcsn_.eps2 = *epsil2;
10153 mmprcsn_.eps3 = *epsil3;
10154 mmprcsn_.eps4 = *epsil4;
10155 mmprcsn_.niterm = *niter1;
10156 mmprcsn_.niterr = *niter2;
10161 //=======================================================================
10162 //function : AdvApp2Var_MathBase::pow__di
10164 //=======================================================================
10165 doublereal AdvApp2Var_MathBase::pow__di (doublereal *x,
10169 register integer ii ;
10170 doublereal result ;
10173 if ( *n > 0 ) {absolute = *n;}
10174 else {absolute = -*n;}
10175 /* System generated locals */
10176 for(ii = 0 ; ii < absolute ; ii++) {
10180 result = 1.0e0 / result ;
10186 /* **********************************************************************
10191 /* Calculate integer function power not obligatory in the most efficient way ;
10198 /* INPUT ARGUMENTS : */
10199 /* ------------------ */
10200 /* X : argument of X**N */
10203 /* OUTPUT ARGUMENTS : */
10204 /* ------------------- */
10207 /* COMMONS USED : */
10208 /* ---------------- */
10210 /* REFERENCES CALLED : */
10211 /* ----------------------- */
10213 /* DESCRIPTION/NOTES/LIMITATIONS : */
10214 /* ----------------------------------- */
10217 /* ***********************************************************************/
10219 //=======================================================================
10220 //function : pow__ii
10222 //=======================================================================
10223 integer pow__ii(integer *x,
10227 register integer ii ;
10231 if ( *n > 0 ) {absolute = *n;}
10232 else {absolute = -*n;}
10233 /* System generated locals */
10234 for(ii = 0 ; ii < absolute ; ii++) {
10238 result = 1 / result ;
10244 /* **********************************************************************
10246 /* **********************************************************************
10251 /* Calculate integer function power not obligatory in the most efficient way ;
10258 /* INPUT ARGUMENTS : */
10259 /* ------------------ */
10260 /* X : argument of X**N */
10263 /* OUTPUT ARGUMENTS : */
10264 /* ------------------- */
10267 /* COMMONS USED : */
10268 /* ---------------- */
10270 /* REFERENCES CALLED : */
10271 /* ----------------------- */
10273 /* DESCRIPTION/NOTES/LIMITATIONS : */
10274 /* ----------------------------------- */
10277 /* ***********************************************************************/
10279 //=======================================================================
10280 //function : AdvApp2Var_MathBase::msc_
10282 //=======================================================================
10283 doublereal AdvApp2Var_MathBase::msc_(integer *ndimen,
10284 doublereal *vecte1,
10285 doublereal *vecte2)
10288 /* System generated locals */
10290 doublereal ret_val;
10292 /* Local variables */
10293 static integer i__;
10294 static doublereal x;
10298 /************************************************************************
10303 /* Calculate the scalar product of 2 vectors in the space */
10304 /* of dimension NDIMEN. */
10308 /* PRODUCT MSCALAIRE. */
10310 /* INPUT ARGUMENTS : */
10311 /* ------------------ */
10312 /* NDIMEN : Dimension of the space. */
10313 /* VECTE1,VECTE2: Vectors. */
10315 /* OUTPUT ARGUMENTS : */
10316 /* ------------------- */
10318 /* COMMONS USED : */
10319 /* ---------------- */
10321 /* REFERENCES CALLED : */
10322 /* ----------------------- */
10324 /* DESCRIPTION/NOTES/LIMITATIONS : */
10325 /* ----------------------------------- */
10328 /* ***********************************************************************
10332 /* PRODUIT MSCALAIRE */
10333 /* Parameter adjustments */
10337 /* Function Body */
10341 for (i__ = 1; i__ <= i__1; ++i__) {
10342 x += vecte1[i__] * vecte2[i__];
10347 /* ----------------------------------- THE END --------------------------
10353 //=======================================================================
10354 //function : mvcvin2_
10356 //=======================================================================
10357 int mvcvin2_(integer *ncoeff,
10358 doublereal *crvold,
10359 doublereal *crvnew,
10363 /* System generated locals */
10364 integer i__1, i__2;
10366 /* Local variables */
10367 static integer m1jm1, ncfm1, j, k;
10368 static doublereal bid;
10369 static doublereal cij1, cij2;
10373 /************************************************************************
10378 /* INVERSION OF THE PARAMETERS ON CURVE 2D. */
10382 /* CURVE,2D,INVERSION,PARAMETER. */
10384 /* INPUT ARGUMENTS : */
10385 /* ------------------ */
10386 /* NCOEFF : NB OF COEFF OF THE CURVE. */
10387 /* CRVOLD : CURVE OF ORIGIN */
10389 /* OUTPUT ARGUMENTS : */
10390 /* ------------------- */
10391 /* CRVNEW : THE RESULTING CURVE AFTER CHANGE OF T BY 1-T */
10392 /* IERCOD : 0 OK, */
10393 /* 10 NB OF COEFF NULL OR TOO GREAT. */
10395 /* COMMONS USED : */
10396 /* ---------------- */
10399 /* REFERENCES CALLED : */
10400 /* ---------------------- */
10402 /* DESCRIPTION/NOTES/LIMITATIONS : */
10403 /* ----------------------------------- */
10404 /* THE FOLLOWING CALL IS ABSOLUTELY LEGAL : */
10405 /* CALL MVCVIN2(NCOEFF,CURVE,CURVE,IERCOD), THE TABLE CURVE */
10406 /* BECOMES INPUT AND OUTPUT ARGUMENT (RBD). */
10407 /* BECAUSE OF MCCNP, THE NB OF COEFF OF THE CURVE IS LIMITED TO */
10408 /* NDGCNP+1 = 61. */
10411 /* ***********************************************************************
10415 /* **********************************************************************
10420 /* Serves to provide coefficients of the binome (triangle of Pascal). */
10424 /* Coeff of binome from 0 to 60. read only . init par block data */
10426 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
10427 /* ----------------------------------- */
10428 /* The coefficients of the binome form a triangular matrix. */
10429 /* This matrix is completed in table CNP by transposition. */
10430 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
10432 /* Initialization is done by block-data MMLLL09.RES, */
10433 /* created by program MQINICNP.FOR (see the team (AC) ). */
10437 /* **********************************************************************
10442 /* ***********************************************************************
10445 /* Parameter adjustments */
10449 /* Function Body */
10450 if (*ncoeff < 1 || *ncoeff - 1 > 60) {
10457 /* CONSTANT TERM OF THE NEW CURVE */
10462 for (k = 2; k <= i__1; ++k) {
10463 cij1 += crvold[(k << 1) + 1];
10464 cij2 += crvold[(k << 1) + 2];
10468 if (*ncoeff == 1) {
10472 /* INTERMEDIARY POWERS OF THE PARAMETER */
10474 ncfm1 = *ncoeff - 1;
10477 for (j = 2; j <= i__1; ++j) {
10479 cij1 = crvold[(j << 1) + 1];
10480 cij2 = crvold[(j << 1) + 2];
10482 for (k = j + 1; k <= i__2; ++k) {
10483 bid = mmcmcnp_.cnp[k - 1 + (j - 1) * 61];
10484 cij1 += crvold[(k << 1) + 1] * bid;
10485 cij2 += crvold[(k << 1) + 2] * bid;
10487 crvnew[(j << 1) + 1] = cij1 * m1jm1;
10488 crvnew[(j << 1) + 2] = cij2 * m1jm1;
10491 /* TERM OF THE HIGHEST DEGREE */
10493 crvnew[(*ncoeff << 1) + 1] = -crvold[(*ncoeff << 1) + 1] * m1jm1;
10494 crvnew[(*ncoeff << 1) + 2] = -crvold[(*ncoeff << 1) + 2] * m1jm1;
10498 AdvApp2Var_SysBase::maermsg_("MVCVIN2", iercod, 7L);
10503 //=======================================================================
10504 //function : mvcvinv_
10506 //=======================================================================
10507 int mvcvinv_(integer *ncoeff,
10508 doublereal *crvold,
10509 doublereal *crvnew,
10513 /* System generated locals */
10514 integer i__1, i__2;
10516 /* Local variables */
10517 static integer m1jm1, ncfm1, j, k;
10518 static doublereal bid;
10519 //extern /* Subroutine */ int maermsg_();
10520 static doublereal cij1, cij2, cij3;
10523 /* **********************************************************************
10528 /* INVERSION OF THE PARAMETER ON A CURBE 3D (I.E. INVERSION */
10529 /* OF THE DIRECTION OF PARSING). */
10533 /* CURVE,INVERSION,PARAMETER. */
10535 /* INPUT ARGUMENTS : */
10536 /* ------------------ */
10537 /* NCOEFF : NB OF COEFF OF THE CURVE. */
10538 /* CRVOLD : CURVE OF ORIGIN */
10540 /* OUTPUT ARGUMENTS : */
10541 /* ------------------- */
10542 /* CRVNEW : RESULTING CURVE AFTER CHANGE OF T INTO 1-T */
10543 /* IERCOD : 0 OK, */
10544 /* 10 NB OF COEFF NULL OR TOO GREAT. */
10546 /* COMMONS USED : */
10547 /* ---------------- */
10550 /* REFERENCES CALLED : */
10551 /* ---------------------- */
10553 /* DESCRIPTION/NOTES/LIMITATIONS : */
10554 /* ----------------------------------- */
10555 /* THE FOLLOWING CALL IS ABSOLUTELY LEGAL : */
10556 /* CALL MVCVINV(NCOEFF,CURVE,CURVE,IERCOD), TABLE CURVE */
10557 /* BECOMES INPUT AND OUTPUT ARGUMENT (RBD). */
10558 /* THE NUMBER OF COEFF OF THE CURVE IS LIMITED TO NDGCNP+1 = 61 */
10559 /* BECAUSE OF USE OF COMMON MCCNP. */
10561 /* ***********************************************************************
10564 /* **********************************************************************
10569 /* Serves to provide the binomial coefficients (triangle of Pascal). */
10573 /* Binomial Coeff from 0 to 60. read only . init par block data */
10575 /* DEMSCRIPTION/NOTES/LIMITATIONS : */
10576 /* ----------------------------------- */
10577 /* The binomial coefficients form a triangular matrix. */
10578 /* This matrix is completed in table CNP by its transposition. */
10579 /* So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
10581 /* Initialisation is done by block-data MMLLL09.RES, */
10582 /* created by program MQINICNP.FOR (see the team (AC) ). */
10584 /* **********************************************************************
10589 /* ***********************************************************************
10592 /* Parameter adjustments */
10596 /* Function Body */
10597 if (*ncoeff < 1 || *ncoeff - 1 > 60) {
10603 /* CONSTANT TERM OF THE NEW CURVE */
10609 for (k = 2; k <= i__1; ++k) {
10610 cij1 += crvold[k * 3 + 1];
10611 cij2 += crvold[k * 3 + 2];
10612 cij3 += crvold[k * 3 + 3];
10618 if (*ncoeff == 1) {
10622 /* INTERMEDIARY POWER OF THE PARAMETER */
10624 ncfm1 = *ncoeff - 1;
10627 for (j = 2; j <= i__1; ++j) {
10629 cij1 = crvold[j * 3 + 1];
10630 cij2 = crvold[j * 3 + 2];
10631 cij3 = crvold[j * 3 + 3];
10633 for (k = j + 1; k <= i__2; ++k) {
10634 bid = mmcmcnp_.cnp[k - 1 + (j - 1) * 61];
10635 cij1 += crvold[k * 3 + 1] * bid;
10636 cij2 += crvold[k * 3 + 2] * bid;
10637 cij3 += crvold[k * 3 + 3] * bid;
10640 crvnew[j * 3 + 1] = cij1 * m1jm1;
10641 crvnew[j * 3 + 2] = cij2 * m1jm1;
10642 crvnew[j * 3 + 3] = cij3 * m1jm1;
10646 /* TERM OF THE HIGHEST DEGREE */
10648 crvnew[*ncoeff * 3 + 1] = -crvold[*ncoeff * 3 + 1] * m1jm1;
10649 crvnew[*ncoeff * 3 + 2] = -crvold[*ncoeff * 3 + 2] * m1jm1;
10650 crvnew[*ncoeff * 3 + 3] = -crvold[*ncoeff * 3 + 3] * m1jm1;
10653 AdvApp2Var_SysBase::maermsg_("MVCVINV", iercod, 7L);
10657 //=======================================================================
10658 //function : mvgaus0_
10660 //=======================================================================
10661 int mvgaus0_(integer *kindic,
10662 doublereal *urootl,
10663 doublereal *hiltab,
10668 /* System generated locals */
10671 /* Local variables */
10672 static doublereal tamp[40];
10673 static integer ndegl, kg, ii;
10675 /* **********************************************************************
10680 /* Loading of a degree gives roots of LEGENDRE polynom */
10681 /* DEFINED on [-1,1] and weights of Gauss quadrature formulas */
10682 /* (based on corresponding LAGRANGIAN interpolators). */
10683 /* The symmetry relative to 0 is used between [-1,0] and [0,1]. */
10687 /* . VOLUMIC, LEGENDRE, LAGRANGE, GAUSS */
10689 /* INPUT ARGUMENTSE : */
10690 /* ------------------ */
10692 /* KINDIC : Takes values from 1 to 10 depending of the degree */
10693 /* of the used polynom. */
10694 /* The degree of the polynom is equal to 4 k, i.e. 4, 8, */
10695 /* 12, 16, 20, 24, 28, 32, 36 and 40. */
10697 /* OUTPUT ARGUMENTS : */
10698 /* ------------------- */
10700 /* UROOTL : Roots of LEGENDRE polynom in domain [1,0] */
10701 /* given in decreasing order. For domain [-1,0], it is */
10702 /* necessary to take the opposite values. */
10703 /* HILTAB : LAGRANGE interpolators associated to roots. For */
10704 /* opposed roots, interpolatorsare equal. */
10705 /* NBRVAL : Nb of coefficients. Is equal to the half of degree */
10706 /* depending on the symmetry (i.e. 2*KINDIC). */
10708 /* IERCOD : Error code: */
10709 /* < 0 ==> Attention - Warning */
10710 /* =-1 ==> Value of false KINDIC. NBRVAL is forced to 20 */
10712 /* = 0 ==> Everything is OK */
10714 /* COMMON USED : */
10715 /* ---------------- */
10717 /* REFERENCES CALLED : */
10718 /* ------------------- */
10720 /* DESCRIPTION/NOTES/LIMITATIONS : */
10721 /* --------------------------------- */
10722 /* If KINDIC is not correct (i.e < 1 or > 10), the degree is set */
10723 /* to 40 directly (ATTENTION to overload - to avoid it, */
10724 /* preview UROOTL and HILTAB dimensioned at least to 20). */
10726 /* The value of coefficients was calculated with quadruple precision
10727 /* by JJM with help of GD. */
10728 /* Checking of roots was done by GD. */
10730 /* See detailed explications on the listing */
10732 /* **********************************************************************
10736 /* ------------------------------------ */
10737 /* ****** Test validity of KINDIC ** */
10738 /* ------------------------------------ */
10740 /* Parameter adjustments */
10744 /* Function Body */
10747 if (kg < 1 || kg > 10) {
10752 ndegl = *nbrval << 1;
10754 /* ----------------------------------------------------------------------
10756 /* ****** Load NBRVAL positive roots depending on the degree **
10758 /* ----------------------------------------------------------------------
10760 /* ATTENTION : Sign minus (-) in the loop is intentional. */
10762 mmextrl_(&ndegl, tamp);
10764 for (ii = 1; ii <= i__1; ++ii) {
10765 urootl[ii] = -tamp[ii - 1];
10769 /* ------------------------------------------------------------------- */
10770 /* ****** Loading of NBRVAL Gauss weight depending on the degree ** */
10771 /* ------------------------------------------------------------------- */
10773 mmexthi_(&ndegl, tamp);
10775 for (ii = 1; ii <= i__1; ++ii) {
10776 hiltab[ii] = tamp[ii - 1];
10780 /* ------------------------------- */
10781 /* ****** End of sub-program ** */
10782 /* ------------------------------- */
10787 //=======================================================================
10788 //function : mvpscr2_
10790 //=======================================================================
10791 int mvpscr2_(integer *ncoeff,
10792 doublereal *curve2,
10793 doublereal *tparam,
10794 doublereal *pntcrb)
10796 /* System generated locals */
10799 /* Local variables */
10800 static integer ndeg, kk;
10801 static doublereal xxx, yyy;
10805 /* **********************************************************************
10810 /* POSITIONING ON CURVE (NCF,2) IN SPACE OF DIMENSION 2. */
10814 /* TOUS,MATH_ACCES:: COURBE&,POSITIONNEMENT,&POINT. */
10816 /* INPUT ARGUMENTS : */
10817 /* ------------------ */
10818 /* NCOEFF : NUMBER OF COEFFICIENTS OF THE CURVE */
10819 /* CURVE2 : EQUATION OF CURVE 2D */
10820 /* TPARAM : VALUE OF PARAMETER AT GIVEN POINT */
10822 /* OUTPUT ARGUMENTS : */
10823 /* ------------------- */
10824 /* PNTCRB : COORDINATES OF POINT CORRESPONDING TO PARAMETER */
10825 /* TPARAM ON CURVE 2D CURVE2. */
10827 /* COMMONS USED : */
10828 /* ---------------- */
10830 /* REFERENCES CALLED : */
10831 /* ---------------------- */
10833 /* DESCRIPTION/NOTES/LIMITATIONS : */
10834 /* ----------------------------------- */
10835 /* MSCHEMA OF HORNER. */
10838 /* **********************************************************************
10842 /* -------- INITIALIZATIONS AND PROCESSING OF PARTICULAR CASES ----------
10845 /* ---> Cas when NCOEFF > 1 (case STANDARD). */
10846 /* Parameter adjustments */
10850 /* Function Body */
10851 if (*ncoeff >= 2) {
10854 /* ---> Case when NCOEFF <= 1. */
10855 if (*ncoeff <= 0) {
10859 } else if (*ncoeff == 1) {
10860 pntcrb[1] = curve2[3];
10861 pntcrb[2] = curve2[4];
10865 /* -------------------- MSCHEMA OF HORNER (PARTICULAR CASE) --------------
10870 if (*tparam == 1.) {
10874 for (kk = 1; kk <= i__1; ++kk) {
10875 xxx += curve2[(kk << 1) + 1];
10876 yyy += curve2[(kk << 1) + 2];
10880 } else if (*tparam == 0.) {
10881 pntcrb[1] = curve2[3];
10882 pntcrb[2] = curve2[4];
10886 /* ---------------------------- MSCHEMA OF HORNER ------------------------
10888 /* ---> TPARAM is different from 1.D0 and 0.D0. */
10890 ndeg = *ncoeff - 1;
10891 xxx = curve2[(*ncoeff << 1) + 1];
10892 yyy = curve2[(*ncoeff << 1) + 2];
10893 for (kk = ndeg; kk >= 1; --kk) {
10894 xxx = xxx * *tparam + curve2[(kk << 1) + 1];
10895 yyy = yyy * *tparam + curve2[(kk << 1) + 2];
10900 /* ------------------------ RECOVER THE CALCULATED POINT ---------------
10907 /* ------------------------------ THE END -------------------------------
10914 //=======================================================================
10915 //function : mvpscr3_
10917 //=======================================================================
10918 int mvpscr3_(integer *ncoeff,
10919 doublereal *curve3,
10920 doublereal *tparam,
10921 doublereal *pntcrb)
10924 /* System generated locals */
10927 /* Local variables */
10928 static integer ndeg, kk;
10929 static doublereal xxx, yyy, zzz;
10933 /* **********************************************************************
10938 /* POSITIONING ON A CURVE (3,NCF) IN THE SPACE OF DIMENSION 3. */
10942 /* TOUS, MATH_ACCES:: COURBE&,POSITIONNEMENT,&POINT. */
10944 /* INPUT ARGUMENTS : */
10945 /* ------------------ */
10946 /* NCOEFF : NB OF COEFFICIENTS OF THE CURVE */
10947 /* CURVE3 : EQUATION OF CURVE 3D */
10948 /* TPARAM : VALUE OF THE PARAMETER AT THE GIVEN POINT */
10950 /* OUTPUT ARGUMENTS : */
10951 /* ------------------- */
10952 /* PNTCRB : COORDINATES OF THE POINT CORRESPONDING TO PARAMETER */
10953 /* TPARAM ON CURVE 3D CURVE3. */
10955 /* COMMONS USED : */
10956 /* ---------------- */
10958 /* REFERENCES CALLED : */
10959 /* ---------------------- */
10962 /* DESCRIPTION/NOTES/LIMITATIONS : */
10963 /* ----------------------------------- */
10964 /* MSCHEMA OF HORNER. */
10966 /* **********************************************************************
10969 /* **********************************************************************
10973 /* -------- INITIALISATIONS AND PROCESSING OF PARTICULAR CASES ----------
10976 /* ---> Case when NCOEFF > 1 (cas STANDARD). */
10977 /* Parameter adjustments */
10981 /* Function Body */
10982 if (*ncoeff >= 2) {
10985 /* ---> Case when NCOEFF <= 1. */
10986 if (*ncoeff <= 0) {
10991 } else if (*ncoeff == 1) {
10992 pntcrb[1] = curve3[4];
10993 pntcrb[2] = curve3[5];
10994 pntcrb[3] = curve3[6];
10998 /* -------------------- MSCHEMA OF HORNER (PARTICULAR CASE) --------------
11003 if (*tparam == 1.) {
11008 for (kk = 1; kk <= i__1; ++kk) {
11009 xxx += curve3[kk * 3 + 1];
11010 yyy += curve3[kk * 3 + 2];
11011 zzz += curve3[kk * 3 + 3];
11015 } else if (*tparam == 0.) {
11016 pntcrb[1] = curve3[4];
11017 pntcrb[2] = curve3[5];
11018 pntcrb[3] = curve3[6];
11022 /* ---------------------------- MSCHEMA OF HORNER ------------------------
11024 /* ---> Here TPARAM is different from 1.D0 and 0.D0. */
11026 ndeg = *ncoeff - 1;
11027 xxx = curve3[*ncoeff * 3 + 1];
11028 yyy = curve3[*ncoeff * 3 + 2];
11029 zzz = curve3[*ncoeff * 3 + 3];
11030 for (kk = ndeg; kk >= 1; --kk) {
11031 xxx = xxx * *tparam + curve3[kk * 3 + 1];
11032 yyy = yyy * *tparam + curve3[kk * 3 + 2];
11033 zzz = zzz * *tparam + curve3[kk * 3 + 3];
11038 /* ------------------------ RETURN THE CALCULATED POINT ------------------
11046 /* ------------------------------ THE END -------------------------------
11053 //=======================================================================
11054 //function : AdvApp2Var_MathBase::mvsheld_
11056 //=======================================================================
11057 int AdvApp2Var_MathBase::mvsheld_(integer *n,
11063 /* System generated locals */
11064 integer dtab_dim1, dtab_offset, i__1, i__2;
11066 /* Local variables */
11067 static integer incr;
11068 static doublereal dsave;
11069 static integer i3, i4, i5, incrp1;
11072 /************************************************************************
11077 /* PARSING OF COLUMNS OF TABLE OF REAL*8 BY SHELL METHOD*/
11078 /* (IN INCREASING ORDER) */
11082 /* POINT-ENTRY, PARSING, SHELL */
11084 /* INPUT ARGUMENTS : */
11085 /* ------------------ */
11086 /* N : NUMBER OF COLUMNS OF THE TABLE */
11087 /* IS : NUMBER OF LINE OF THE TABLE */
11088 /* DTAB : TABLE OF REAL*8 TO BE PARSED */
11089 /* ICLE : POSITION OF THE KEY ON THE COLUMN */
11091 /* OUTPUT ARGUMENTS : */
11092 /* ------------------- */
11093 /* DTAB : PARSED TABLE */
11095 /* COMMONS USED : */
11096 /* ---------------- */
11099 /* REFERENCES CALLED : */
11100 /* ---------------------- */
11103 /* DESCRIPTION/NOTES/LIMITATIONS : */
11104 /* ----------------------------------- */
11105 /* CLASSIC SHELL METHOD : PARSING BY SERIES */
11106 /* Declaration DTAB(IS, 1) corresponds to DTAB(IS, *) */
11108 /* ***********************************************************************
11112 /* Parameter adjustments */
11114 dtab_offset = dtab_dim1 + 1;
11115 dtab -= dtab_offset;
11117 /* Function Body */
11121 /* ------------------------ */
11123 /* INITIALIZATION OF THE SEQUENCE OF INCREMENTS */
11124 /* FIND THE GREATEST INCREMENT SO THAT INCR < N/9 */
11128 if (incr >= *n / 9) {
11131 /* ----------------------------- */
11132 incr = incr * 3 + 1;
11135 /* LOOP ON INCREMENTS TILL INCR = 1 */
11136 /* PARSING BY SERIES DISTANT FROM INCR */
11140 /* ----------------- */
11142 for (i3 = incrp1; i3 <= i__1; ++i3) {
11143 /* ---------------------- */
11145 /* SET ELEMENT I3 AT ITS PLACE IN THE SERIES */
11152 /* ------------------------- */
11153 if (dtab[*icle + i4 * dtab_dim1] <= dtab[*icle + (i4 + incr) *
11159 for (i5 = 1; i5 <= i__2; ++i5) {
11160 /* ------------------ */
11161 dsave = dtab[i5 + i4 * dtab_dim1];
11162 dtab[i5 + i4 * dtab_dim1] = dtab[i5 + (i4 + incr) * dtab_dim1];
11163 dtab[i5 + (i4 + incr) * dtab_dim1] = dsave;
11174 /* PASSAGE TO THE NEXT INCREMENT */
11185 //=======================================================================
11186 //function : AdvApp2Var_MathBase::mzsnorm_
11188 //=======================================================================
11189 doublereal AdvApp2Var_MathBase::mzsnorm_(integer *ndimen,
11190 doublereal *vecteu)
11193 /* System generated locals */
11195 doublereal ret_val, d__1, d__2;
11197 /* Local variables */
11198 static doublereal xsom;
11199 static integer i__, irmax;
11203 /* ***********************************************************************
11208 /* SERVES to calculate the euclidian norm of a vector : */
11209 /* ____________________________ */
11210 /* Z = V V(1)**2 + V(2)**2 + ... */
11216 /* INPUT ARGUMENTS : */
11217 /* ------------------ */
11218 /* NDIMEN : Dimension of the vector */
11219 /* VECTEU : vector of dimension NDIMEN */
11221 /* OUTPUT ARGUMENTS : */
11222 /* ------------------- */
11223 /* MZSNORM : Value of the euclidian norm of vector VECTEU */
11225 /* COMMONS USED : */
11226 /* ---------------- */
11230 /* REFERENCES CALLED : */
11231 /* ---------------------- */
11233 /* R*8 ABS R*8 SQRT */
11235 /* DESCRIPTION/NOTESS/LIMITATIONS : */
11236 /* ----------------------------------- */
11237 /* To limit the risks of overflow, */
11238 /* the term of the strongest absolute value is factorized : */
11239 /* _______________________ */
11240 /* Z = !V(1)! * V 1 + (V(2)/V(1))**2 + ... */
11243 /* ***********************************************************************
11246 /* ***********************************************************************
11250 /* ***********************************************************************
11253 /* ***********************************************************************
11256 /* ___ Find the strongest absolute value term */
11258 /* Parameter adjustments */
11261 /* Function Body */
11264 for (i__ = 2; i__ <= i__1; ++i__) {
11265 if ((d__1 = vecteu[irmax], advapp_abs(d__1)) < (d__2 = vecteu[i__], advapp_abs(d__2)
11272 /* ___ Calculate the norme */
11274 if ((d__1 = vecteu[irmax], advapp_abs(d__1)) < 1.) {
11277 for (i__ = 1; i__ <= i__1; ++i__) {
11278 /* Computing 2nd power */
11279 d__1 = vecteu[i__];
11280 xsom += d__1 * d__1;
11283 ret_val = sqrt(xsom);
11287 for (i__ = 1; i__ <= i__1; ++i__) {
11288 if (i__ == irmax) {
11291 /* Computing 2nd power */
11292 d__1 = vecteu[i__] / vecteu[irmax];
11293 xsom += d__1 * d__1;
11297 ret_val = (d__1 = vecteu[irmax], advapp_abs(d__1)) * sqrt(xsom);
11300 /* ***********************************************************************
11302 /* RETURN CALLING PROGRAM */
11303 /* ***********************************************************************