2 // AdvApp2Var_ApproxF2var.cxx
5 #include <AdvApp2Var_SysBase.hxx>
6 #include <AdvApp2Var_MathBase.hxx>
7 #include <AdvApp2Var_Data_f2c.hxx>
8 #include <AdvApp2Var_Data.hxx>
9 #include <AdvApp2Var_ApproxF2var.hxx>
13 int mmjacpt_(const integer *ndimen,
14 const integer *ncoefu,
15 const integer *ncoefv,
16 const integer *iordru,
17 const integer *iordrv,
18 const doublereal *ptclgd,
25 int mma2ce2_(integer *numdec,
60 int mma2cfu_(integer *ndujac,
72 int mma2cfv_(integer *ndvjac,
82 int mma2er1_(integer *ndjacu,
98 int mma2er2_(integer *ndjacu,
117 int mma2moy_(integer *ndgumx,
130 int mma2ds2_(integer *ndimen,
163 int mma1fdi_(integer *ndimen,
165 void (*foncnp) (// see AdvApp2Var_EvaluatorFunc2Var.hxx for details
192 int mma1cdi_(integer *ndimen,
204 int mma1jak_(integer *ndimen,
214 int mma1cnt_(integer *ndimen,
223 int mma1fer_(integer *ndimen,
238 int mma1noc_(doublereal *dfuvin,
249 int mmmapcoe_(integer *ndim,
259 int mmaperm_(integer *ncofmx,
268 #define mmapgss_1 mmapgss_
269 #define mmapgs0_1 mmapgs0_
270 #define mmapgs1_1 mmapgs1_
271 #define mmapgs2_1 mmapgs2_
273 //=======================================================================
274 //function : mma1cdi_
276 //=======================================================================
277 int mma1cdi_(integer *ndimen,
289 static integer c__1 = 1;
291 /* System generated locals */
292 integer contr1_dim1, contr1_offset, contr2_dim1, contr2_offset,
293 somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
294 fpntab_dim1, fpntab_offset, hermit_dim1, hermit_offset, i__1,
297 /* Local variables */
298 static integer nroo2, ncfhe, nd, ii, kk;
299 static integer ibb, kkm, kkp;
300 static doublereal bid1, bid2, bid3;
302 /* **********************************************************************
306 /* Discretisation on the parameters of interpolation polynomes */
307 /* constraints of order IORDRE. */
311 /* ALL, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
313 /* INPUT ARGUMENTS : */
314 /* ------------------ */
315 /* NDIMEN: Space dimension. */
316 /* NBROOT: Number of INTERNAL discretisation parameters. */
317 /* It is also the root number Legendre polynome where */
318 /* the discretization is performed. */
319 /* ROOTLG: Table of discretization parameters ON (-1,1). */
320 /* IORDRE: Order of constraint imposed to the extremities of the iso. */
321 /* = 0, the extremities of the iso are calculated */
322 /* = 1, additionally, the 1st derivative in the direction */
323 /* of the iso is calculated. */
324 /* = 2, additionally, the 2nd derivative in the direction */
325 /* of the iso is calculated. */
326 /* CONTR1: Contains, if IORDRE>=0, values IORDRE+1 in TTABLE(0)
328 /* (1st extremity) of derivatives of F(Uc,Ve) or F(Ue,Vc), */
330 /* CONTR2: Contains, if IORDRE>=0, values IORDRE+1 in */
331 /* TTABLE(NBROOT+1) (2nd extremity) of: */
332 /* If ISOFAV=1, derived of order IDERIV by U, derived */
333 /* ordre 0 to IORDRE by V of F(Uc,Ve) or Uc=TCONST */
334 /* (fixed iso value) and Ve is the fixed extremity. */
335 /* If ISOFAV=2, derivative of order IDERIV by V, derivative */
336 /* of order 0 to IORDRE by U of F(Ue,Vc) or Vc=TCONST */
337 /* (fixed iso value) and Ue is the fixed extremity. */
339 /* SOMTAB: Table of NBROOT/2 sums of 2 index points */
340 /* NBROOT-II+1 and II, for II = 1, NBROOT/2. */
341 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
342 /* NBROOT-II+1 and II, for II = 1, NBROOT/2. */
344 /* OUTPUT ARGUMENTS : */
345 /* ------------------- */
346 /* SOMTAB: Table of NBROOT/2 sums of 2 index points */
347 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
348 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
349 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
350 /* FPNTAB: Auxiliary table. */
351 /* HERMIT: Table of coeff. 2*(IORDRE+1) Hermite polynoms */
352 /* of degree 2*IORDRE+1. */
353 /* IERCOD: Error code, */
354 /* = 0, Everythig is OK */
355 /* = 1, The value of IORDRE is out of (0,2) */
357 /* ---------------- */
359 /* REFERENCES CALLED : */
360 /* ----------------------- */
362 /* DESCRIPTION/NOTES/LIMITATIONS : */
363 /* ----------------------------------- */
364 /* The results of discretization are arranged in 2 tables */
365 /* SOMTAB and DIFTAB to earn time during the */
366 /* calculation of coefficients of the approximation curve. */
368 /* If NBROOT is uneven in SOMTAB(0,*) and DIFTAB(0,*) one stores */
369 /* the values of the median root of Legendre (0.D0 in (-1,1)). */
371 /* **********************************************************************
374 /* Name of the routine */
377 /* Parameter adjustments */
378 diftab_dim1 = *nbroot / 2 + 1;
379 diftab_offset = diftab_dim1;
380 diftab -= diftab_offset;
381 somtab_dim1 = *nbroot / 2 + 1;
382 somtab_offset = somtab_dim1;
383 somtab -= somtab_offset;
385 hermit_dim1 = (*iordre << 1) + 2;
386 hermit_offset = hermit_dim1;
387 hermit -= hermit_offset;
388 fpntab_dim1 = *nbroot;
389 fpntab_offset = fpntab_dim1 + 1;
390 fpntab -= fpntab_offset;
391 contr2_dim1 = *ndimen;
392 contr2_offset = contr2_dim1 + 1;
393 contr2 -= contr2_offset;
394 contr1_dim1 = *ndimen;
395 contr1_offset = contr1_dim1 + 1;
396 contr1 -= contr1_offset;
399 ibb = AdvApp2Var_SysBase::mnfndeb_();
401 AdvApp2Var_SysBase::mgenmsg_("MMA1CDI", 7L);
405 /* --- Recuperate 2*(IORDRE+1) coeff of 2*(IORDRE+1) of Hermite polynom ---
408 AdvApp2Var_ApproxF2var::mma1her_(iordre, &hermit[hermit_offset], iercod);
413 /* ------------------- Discretization of Hermite polynoms -----------
416 ncfhe = (*iordre + 1) << 1;
418 for (ii = 1; ii <= i__1; ++ii) {
420 for (kk = 1; kk <= i__2; ++kk) {
421 AdvApp2Var_MathBase::mmmpocur_(&ncfhe, &c__1, &ncfhe, &hermit[ii * hermit_dim1], &
422 rootlg[kk], &fpntab[kk + ii * fpntab_dim1]);
428 /* ---- Discretizations of boundary polynoms are taken ----
433 for (nd = 1; nd <= i__1; ++nd) {
435 for (ii = 1; ii <= i__2; ++ii) {
436 bid1 = contr1[nd + ii * contr1_dim1];
437 bid2 = contr2[nd + ii * contr2_dim1];
439 for (kk = 1; kk <= i__3; ++kk) {
440 kkm = nroo2 - kk + 1;
441 bid3 = bid1 * fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1] +
442 bid2 * fpntab[kkm + (ii << 1) * fpntab_dim1];
443 somtab[kk + nd * somtab_dim1] -= bid3;
444 diftab[kk + nd * diftab_dim1] += bid3;
448 for (kk = 1; kk <= i__3; ++kk) {
449 kkp = (*nbroot + 1) / 2 + kk;
450 bid3 = bid1 * fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] +
451 bid2 * fpntab[kkp + (ii << 1) * fpntab_dim1];
452 somtab[kk + nd * somtab_dim1] -= bid3;
453 diftab[kk + nd * diftab_dim1] -= bid3;
461 /* ------------ Cas when discretization is done on the roots of a -----------
463 /* ---------- Legendre polynom of uneven degree, 0 is root --------
466 if (*nbroot % 2 == 1) {
468 for (nd = 1; nd <= i__1; ++nd) {
470 for (ii = 1; ii <= i__2; ++ii) {
471 bid3 = fpntab[nroo2 + 1 + ((ii << 1) - 1) * fpntab_dim1] *
472 contr1[nd + ii * contr1_dim1] + fpntab[nroo2 + 1 + (
473 ii << 1) * fpntab_dim1] * contr2[nd + ii *
477 somtab[nd * somtab_dim1] -= bid3;
478 diftab[nd * diftab_dim1] -= bid3;
485 /* ------------------------------ The End -------------------------------
487 /* --> IORDRE is not in the authorized zone. */
494 AdvApp2Var_SysBase::mgsomsg_("MMA1CDI", 7L);
499 //=======================================================================
500 //function : mma1cnt_
502 //=======================================================================
503 int mma1cnt_(integer *ndimen,
511 /* System generated locals */
512 integer contr1_dim1, contr1_offset, contr2_dim1, contr2_offset,
513 hermit_dim1, hermit_offset, crvjac_dim1, crvjac_offset, i__1,
516 /* Local variables */
517 static integer nd, ii, jj, ibb;
518 static doublereal bid;
521 /* ***********************************************************************
526 /* Add constraint to polynom. */
530 /* ALL,AB_SPECIFI::COURE&,APPROXIMATION,ADDITION,&CONSTRAINT */
532 /* INPUT ARGUMENTS : */
533 /* -------------------- */
534 /* NDIMEN: Dimension of the space */
535 /* IORDRE: Order of constraint. */
536 /* CONTR1: pt of constraint in -1, from order 0 to IORDRE. */
537 /* CONTR2: Pt of constraint in +1, from order 0 to IORDRE. */
538 /* HERMIT: Table of Hermit polynoms of order IORDRE. */
539 /* CRVJAV: Curve of approximation in Jacobi base. */
541 /* OUTPUT ARGUMENTS : */
542 /* --------------------- */
543 /* CRVJAV: Curve of approximation in Jacobi base */
544 /* to which the polynom of interpolation of constraints is added. */
547 /* ------------------ */
550 /* REFERENCES CALLED : */
551 /* --------------------- */
554 /* DESCRIPTION/NOTES/LIMITATIONS : */
555 /* ----------------------------------- */
558 /* ***********************************************************************
561 /* ***********************************************************************
563 /* Name of the routine */
565 /* ***********************************************************************
567 /* INITIALISATIONS */
568 /* ***********************************************************************
571 /* Parameter adjustments */
572 hermit_dim1 = (*iordre << 1) + 2;
573 hermit_offset = hermit_dim1;
574 hermit -= hermit_offset;
575 contr2_dim1 = *ndimen;
576 contr2_offset = contr2_dim1 + 1;
577 contr2 -= contr2_offset;
578 contr1_dim1 = *ndimen;
579 contr1_offset = contr1_dim1 + 1;
580 contr1 -= contr1_offset;
581 crvjac_dim1 = *ndgjac + 1;
582 crvjac_offset = crvjac_dim1;
583 crvjac -= crvjac_offset;
586 ibb = AdvApp2Var_SysBase::mnfndeb_();
588 AdvApp2Var_SysBase::mgenmsg_("MMA1CNT", 7L);
591 /* ***********************************************************************
594 /* ***********************************************************************
598 for (nd = 1; nd <= i__1; ++nd) {
599 i__2 = (*iordre << 1) + 1;
600 for (ii = 0; ii <= i__2; ++ii) {
603 for (jj = 1; jj <= i__3; ++jj) {
604 bid = bid + contr1[nd + jj * contr1_dim1] *
605 hermit[ii + ((jj << 1) - 1) * hermit_dim1] +
606 contr2[nd + jj * contr2_dim1] * hermit[ii + (jj << 1) * hermit_dim1];
609 crvjac[ii + nd * crvjac_dim1] = bid;
615 /* ***********************************************************************
617 /* RETURN CALLING PROGRAM */
618 /* ***********************************************************************
622 AdvApp2Var_SysBase::mgsomsg_("MMA1CNT", 7L);
628 //=======================================================================
629 //function : mma1fdi_
631 //=======================================================================
632 int mma1fdi_(integer *ndimen,
634 void (*foncnp) (// see AdvApp2Var_EvaluatorFunc2Var.hxx for details
660 /* System generated locals */
661 integer fpntab_dim1, somtab_dim1, somtab_offset, diftab_dim1,
662 diftab_offset, contr1_dim1, contr1_offset, contr2_dim1,
663 contr2_offset, i__1, i__2;
666 /* Local variables */
667 static integer ideb, ifin, nroo2, ideru, iderv;
668 static doublereal renor;
669 static integer ii, nd, ibb, iim, nbp, iip;
670 static doublereal bid1, bid2;
672 /* **********************************************************************
677 /* DiscretiZation of a non-polynomial function F(U,V) or of */
678 /* its derivative with fixed isoparameter. */
682 /* ALL, AB_SPECIFI::FONCTION&, DISCRETISATION, &POINT */
684 /* INPUT ARGUMENTS : */
685 /* ------------------ */
686 /* NDIMEN: Space dimension. */
687 /* UVFONC: Limits of the path of definition by U and by V of the approximated function */
688 /* FONCNP: The NAME of the non-polynomial function to be approximated */
689 /* (external program). */
690 /* ISOFAV: Fixed isoparameter for the discretization; */
691 /* = 1, discretization with fixed U and variable V. */
692 /* = 2, discretization with fixed V and variable U. */
693 /* TCONST: Iso value is also fixed. */
694 /* NBROOT: Number of INTERNAL discretization parameters. */
695 /* (if there are constraints, 2 extremities should be added).
697 /* This is also the root number of the Legendre polynom where */
698 /* the discretization is done. */
699 /* TTABLE: Table of discretization parameters and of 2 extremities */
700 /* (Respectively (-1, NBROOT Legendre roots,1) */
701 /* reframed within the adequate interval. */
702 /* IORDRE: Order of constraint imposed on the extremities of the iso. */
703 /* (If Iso-U, it is necessary to calculate the derivatives by V and vice */
705 /* = 0, the extremities of the iso are calculated. */
706 /* = 1, additionally the 1st derivative in the direction of the iso is calculated */
707 /* = 2, additionally the 2nd derivative in the direction of the iso is calculated */
708 /* IDERIV: Order of derivative transversal to fixed iso (If Iso-U=Uc */
709 /* is fixed, the derivative of order IDERIV is discretized by U of */
710 /* F(Uc,v). Same if iso-V is fixed). */
711 /* Varies from 0 (positioning) to 2 (2nd derivative). */
713 /* OUTPUT ARGUMENTS : */
714 /* ------------------- */
715 /* FPNTAB: Auxiliary table.
716 SOMTAB: Table of NBROOT/2 sums of 2 index points */
717 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
718 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
719 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
720 /* CONTR1: Contains, if IORDRE>=0, values IORDRE+1 in TTABLE(0)
722 /* (1st extremity) of derivatives of F(Uc,Ve) or F(Ue,Vc), */
724 /* CONTR2: Contains, if IORDRE>=0, values IORDRE+1 in */
725 /* TTABLE(NBROOT+1) (2nd extremity) of: */
726 /* If ISOFAV=1, derived of order IDERIV by U, derived */
727 /* ordre 0 to IORDRE by V of F(Uc,Ve) or Uc=TCONST */
728 /* (fixed iso value) and Ve is the fixed extremity. */
729 /* If ISOFAV=2, derivative of order IDERIV by V, derivative */
730 /* of order 0 to IORDRE by U of F(Ue,Vc) or Vc=TCONST */
731 /* (fixed iso value) and Ue is the fixed extremity. */
732 /* IERCOD: Error code > 100; Pb in evaluation of FONCNP, */
733 /* the returned error code is equal to error code of FONCNP + 100. */
736 /* ---------------- */
738 /* REFERENCES CALLED : */
739 /* ----------------------- */
741 /* DESCRIPTION/NOTES/LIMITATIONS : */
742 /* ----------------------------------- */
743 /* The results of discretization are arranged in 2 tables */
744 /* SOMTAB and DIFTAB to earn time during the */
745 /* calculation of coefficients of the approximation curve. */
747 /* If NBROOT is uneven in SOMTAB(0,*) and DIFTAB(0,*) one stores */
748 /* the values of the median root of Legendre (0.D0 in (-1,1)). */
750 /* Function F(u,v) defined in UVFONC is reparameterized in */
751 /* (-1,1)x(-1,1). Then 1st and 2nd derivatives are renormalized. */
754 /* **********************************************************************
757 /* Name of the routine */
760 /* Parameter adjustments */
762 diftab_dim1 = *nbroot / 2 + 1;
763 diftab_offset = diftab_dim1;
764 diftab -= diftab_offset;
765 somtab_dim1 = *nbroot / 2 + 1;
766 somtab_offset = somtab_dim1;
767 somtab -= somtab_offset;
768 fpntab_dim1 = *ndimen;
770 contr2_dim1 = *ndimen;
771 contr2_offset = contr2_dim1 + 1;
772 contr2 -= contr2_offset;
773 contr1_dim1 = *ndimen;
774 contr1_offset = contr1_dim1 + 1;
775 contr1 -= contr1_offset;
778 ibb = AdvApp2Var_SysBase::mnfndeb_();
780 AdvApp2Var_SysBase::mgenmsg_("MMA1FDI", 7L);
784 /* --------------- Definition of the nb of points to calculate --------------
786 /* --> If constraints, the limits are also taken */
790 /* --> Otherwise, only Legendre roots (reframed) are used
796 /* --> Nb of point to calculate. */
797 nbp = ifin - ideb + 1;
800 /* --------------- Determination of the order of global derivation --------
802 /* --> ISOFAV takes only values 1 or 2. */
803 /* if Iso-U, derive by U of order IDERIV */
807 d__1 = (uvfonc[4] - uvfonc[3]) / 2.;
808 renor = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
809 /* if Iso-V, derive by V of order IDERIV */
813 d__1 = (uvfonc[6] - uvfonc[5]) / 2.;
814 renor = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
817 /* ----------- Discretization on roots of the ---------------
819 /* ---------------------- Legendre polynom of degree NBROOT -------------------
831 &fpntab[ideb * fpntab_dim1 + 1],
837 for (nd = 1; nd <= i__1; ++nd) {
839 for (ii = 1; ii <= i__2; ++ii) {
840 iip = (*nbroot + 1) / 2 + ii;
841 iim = nroo2 - ii + 1;
842 bid1 = fpntab[nd + iim * fpntab_dim1];
843 bid2 = fpntab[nd + iip * fpntab_dim1];
844 somtab[ii + nd * somtab_dim1] = renor * (bid2 + bid1);
845 diftab[ii + nd * diftab_dim1] = renor * (bid2 - bid1);
851 /* ------------ Case when discretisation is done on roots of a ----
853 /* ---------- Legendre polynom of uneven degree, 0 is root --------
856 if (*nbroot % 2 == 1) {
858 for (nd = 1; nd <= i__1; ++nd) {
859 somtab[nd * somtab_dim1] = renor * fpntab[nd + (nroo2 + 1) *
861 diftab[nd * diftab_dim1] = renor * fpntab[nd + (nroo2 + 1) *
867 for (nd = 1; nd <= i__1; ++nd) {
868 somtab[nd * somtab_dim1] = 0.;
869 diftab[nd * diftab_dim1] = 0.;
874 /* --------------------- Take into account constraints ----------------
878 /* --> Recover already calculated extremities. */
880 for (nd = 1; nd <= i__1; ++nd) {
881 contr1[nd + contr1_dim1] = renor * fpntab[nd];
882 contr2[nd + contr2_dim1] = renor * fpntab[nd + (*nbroot + 1) *
886 /* --> Nb of points to calculate/call to FONCNP */
888 /* If Iso-U, derive by V till order IORDRE */
890 /* --> Factor of normalisation 1st derivative. */
891 bid1 = (uvfonc[6] - uvfonc[5]) / 2.;
893 for (iderv = 1; iderv <= i__1; ++iderv) {
894 (*foncnp)(ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
895 nbp, ttable, &ideru, &iderv, &contr1[(iderv + 1) *
896 contr1_dim1 + 1], iercod);
903 for (iderv = 1; iderv <= i__1; ++iderv) {
904 (*foncnp)(ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
905 nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
906 iderv + 1) * contr2_dim1 + 1], iercod);
912 /* If Iso-V, derive by U till order IORDRE */
914 /* --> Factor of normalization 1st derivative. */
915 bid1 = (uvfonc[4] - uvfonc[3]) / 2.;
917 for (ideru = 1; ideru <= i__1; ++ideru) {
918 (*foncnp)(ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
919 nbp, ttable, &ideru, &iderv, &contr1[(ideru + 1) *
920 contr1_dim1 + 1], iercod);
927 for (ideru = 1; ideru <= i__1; ++ideru) {
928 (*foncnp)(ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
929 nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
930 ideru + 1) * contr2_dim1 + 1], iercod);
938 /* ------------------------- Normalization of derivatives -------------
940 /* (The function is redefined on (-1,1)*(-1,1)) */
943 for (ii = 1; ii <= i__1; ++ii) {
946 for (nd = 1; nd <= i__2; ++nd) {
947 contr1[nd + (ii + 1) * contr1_dim1] *= bid2;
948 contr2[nd + (ii + 1) * contr2_dim1] *= bid2;
955 /* ------------------------------ The end -------------------------------
961 AdvApp2Var_SysBase::maermsg_("MMA1FDI", iercod, 7L);
964 AdvApp2Var_SysBase::mgsomsg_("MMA1FDI", 7L);
969 //=======================================================================
970 //function : mma1fer_
972 //=======================================================================
973 int mma1fer_(integer *,//ndimen,
987 /* System generated locals */
988 integer crvjac_dim1, crvjac_offset, i__1, i__2;
990 /* Local variables */
991 static integer idim, ncfja, ncfnw, ndses, ii, kk, ibb, ier;
995 /* ***********************************************************************
1000 /* Calculate the degree and the errors of approximation of a border. */
1004 /* TOUS,AB_SPECIFI :: COURBE&,TRONCATURE, &PRECISION */
1006 /* INPUT ARGUMENTS : */
1007 /* -------------------- */
1009 /* NDIMEN: Total Dimension of the space (sum of dimensions of sub-spaces) */
1010 /* NBSESP: Number of "independent" sub-spaces. */
1011 /* NDIMSE: Table of dimensions of sub-spaces. */
1012 /* IORDRE: Order of constraint at the extremities of the border */
1013 /* -1 = no constraints, */
1014 /* 0 = constraints of passage to limits (i.e. C0), */
1015 /* 1 = C0 + constraintes of 1st derivatives (i.e. C1), */
1016 /* 2 = C1 + constraintes of 2nd derivatives (i.e. C2). */
1017 /* NDGJAC: Degree of development in series to use for the calculation
1018 /* in the base of Jacobi. */
1019 /* CRVJAC: Table of coeff. of the curve of approximation in the */
1020 /* base of Jacobi. */
1021 /* NCFLIM: Max number of coeff of the polynomial curve */
1022 /* of approximation (should be above or equal to */
1023 /* 2*IORDRE+2 and below or equal to 50). */
1024 /* EPSAPR: Table of errors of approximations that cannot be passed, */
1025 /* sub-space by sub-space. */
1027 /* OUTPUT ARGUMENTS : */
1028 /* --------------------- */
1029 /* YCVMAX: Auxiliary Table. */
1030 /* ERRMAX: Table of errors (sub-space by sub-space) */
1031 /* MAXIMUM made in the approximation of FONCNP by */
1033 /* ERRMOY: Table of errors (sub-space by sub-space) */
1034 /* AVERAGE made in the approximation of FONCNP by */
1036 /* NCOEFF: Number of significative coeffs. of the calculated "curve". */
1037 /* IERCOD: Error code */
1039 /* =-1, warning, required tolerance can't be */
1040 /* met with coefficients NFCLIM. */
1041 /* = 1, order of constraints (IORDRE) is not within authorised values */
1044 /* COMMONS USED : */
1045 /* ------------------ */
1047 /* REFERENCES CALLED : */
1048 /* --------------------- */
1050 /* DESCRIPTION/NOTES/LIMITATIONS : */
1051 /* ----------------------------------- */
1053 /* **********************************************************************
1056 /* Name of the routine */
1059 /* Parameter adjustments */
1065 crvjac_dim1 = *ndgjac + 1;
1066 crvjac_offset = crvjac_dim1;
1067 crvjac -= crvjac_offset;
1070 ibb = AdvApp2Var_SysBase::mnfndeb_();
1072 AdvApp2Var_SysBase::mgenmsg_("MMA1FER", 7L);
1077 ncfja = *ndgjac + 1;
1079 /* ------------ Calculate the degree of the curve and of the Max error --------
1081 /* -------------- of approximation for all sub-spaces --------
1085 for (ii = 1; ii <= i__1; ++ii) {
1088 /* ------------ cutting of coeff. and calculation of Max error -------
1091 AdvApp2Var_MathBase::mmtrpjj_(&ncfja, &ndses, &ncfja, &epsapr[ii], iordre, &crvjac[idim *
1092 crvjac_dim1], &ycvmax[1], &errmax[ii], &ncfnw);
1094 /* ******************************************************************
1096 /* ------------- If precision OK, calculate the average error -------
1098 /* ******************************************************************
1101 if (ncfnw <= *ncflim) {
1102 mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1103 crvjac_dim1], &ncfnw, &errmoy[ii]);
1104 *ncoeff = max(ncfnw,*ncoeff);
1106 /* ------------- Set the declined coefficients to 0.D0 -----------
1109 nbr0 = *ncflim - ncfnw;
1112 for (kk = 1; kk <= i__2; ++kk) {
1113 AdvApp2Var_SysBase::mvriraz_(&nbr0,
1114 (char *)&crvjac[ncfnw + (idim + kk - 1) * crvjac_dim1]);
1120 /* **************************************************************
1122 /* ------------------- If required precision can't be reached----
1124 /* **************************************************************
1129 /* ------------------------- calculate the Max error ------------
1132 AdvApp2Var_MathBase::mmaperx_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1133 crvjac_dim1], ncflim, &ycvmax[1], &errmax[ii], &ier);
1138 /* -------------------- nb of coeff to be returned -------------
1143 /* ------------------- and calculation of the average error ----
1146 mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1147 crvjac_dim1], ncflim, &errmoy[ii]);
1155 /* ------------------------------ The end -------------------------------
1157 /* --> The order of constraints is not within autorized values. */
1164 AdvApp2Var_SysBase::maermsg_("MMA1FER", iercod, 7L);
1167 AdvApp2Var_SysBase::mgsomsg_("MMA1FER", 7L);
1173 //=======================================================================
1174 //function : mma1her_
1176 //=======================================================================
1177 int AdvApp2Var_ApproxF2var::mma1her_(const integer *iordre,
1181 /* System generated locals */
1182 integer hermit_dim1, hermit_offset;
1184 /* Local variables */
1189 /* **********************************************************************
1194 /* Calculate 2*(IORDRE+1) Hermit polynoms of degree 2*IORDRE+1 */
1199 /* ALL, AB_SPECIFI::CONTRAINTE&, INTERPOLATION, &POLYNOME */
1201 /* INPUT ARGUMENTS : */
1202 /* ------------------ */
1203 /* IORDRE: Order of constraint. */
1204 /* = 0, Polynom of interpolation of order C0 on (-1,1). */
1205 /* = 1, Polynom of interpolation of order C0 and C1 on (-1,1). */
1206 /* = 2, Polynom of interpolation of order C0, C1 and C2 on (-1,1).
1209 /* OUTPUT ARGUMENTS : */
1210 /* ------------------- */
1211 /* HERMIT: Table of 2*IORDRE+2 coeff. of each of 2*(IORDRE+1) */
1212 /* HERMIT polynom. */
1213 /* IERCOD: Error code, */
1215 /* = 1, required order of constraint is not managed here. */
1216 /* COMMONS USED : */
1217 /* ---------------- */
1219 /* REFERENCES CALLED : */
1220 /* ----------------------- */
1222 /* DESCRIPTION/NOTES/LIMITATIONS : */
1223 /* ----------------------------------- */
1224 /* The part of HERMIT(*,2*i+j) table where j=1 or 2 and i=0 to IORDRE,
1225 /* contains the coefficients of the polynom of degree 2*IORDRE+1 */
1226 /* such as ALL values in -1 and in +1 of this polynom and its */
1227 /* derivatives till order of derivation IORDRE are NULL, */
1228 /* EXCEPT for the derivative of order i: */
1229 /* - valued 1 in -1 if j=1 */
1230 /* - valued 1 in +1 if j=2. */
1232 /* **********************************************************************
1235 /* Name of the routine */
1238 /* Parameter adjustments */
1239 hermit_dim1 = (*iordre + 1) << 1;
1240 hermit_offset = hermit_dim1 + 1;
1241 hermit -= hermit_offset;
1244 ibb = AdvApp2Var_SysBase::mnfndeb_();
1246 AdvApp2Var_SysBase::mgenmsg_("MMA1HER", 7L);
1250 /* --- Recover (IORDRE+2) coeff of 2*(IORDRE+1) Hermit polynoms --
1254 hermit[hermit_dim1 + 1] = .5;
1255 hermit[hermit_dim1 + 2] = -.5;
1257 hermit[(hermit_dim1 << 1) + 1] = .5;
1258 hermit[(hermit_dim1 << 1) + 2] = .5;
1259 } else if (*iordre == 1) {
1260 hermit[hermit_dim1 + 1] = .5;
1261 hermit[hermit_dim1 + 2] = -.75;
1262 hermit[hermit_dim1 + 3] = 0.;
1263 hermit[hermit_dim1 + 4] = .25;
1265 hermit[(hermit_dim1 << 1) + 1] = .5;
1266 hermit[(hermit_dim1 << 1) + 2] = .75;
1267 hermit[(hermit_dim1 << 1) + 3] = 0.;
1268 hermit[(hermit_dim1 << 1) + 4] = -.25;
1270 hermit[hermit_dim1 * 3 + 1] = .25;
1271 hermit[hermit_dim1 * 3 + 2] = -.25;
1272 hermit[hermit_dim1 * 3 + 3] = -.25;
1273 hermit[hermit_dim1 * 3 + 4] = .25;
1275 hermit[(hermit_dim1 << 2) + 1] = -.25;
1276 hermit[(hermit_dim1 << 2) + 2] = -.25;
1277 hermit[(hermit_dim1 << 2) + 3] = .25;
1278 hermit[(hermit_dim1 << 2) + 4] = .25;
1279 } else if (*iordre == 2) {
1280 hermit[hermit_dim1 + 1] = .5;
1281 hermit[hermit_dim1 + 2] = -.9375;
1282 hermit[hermit_dim1 + 3] = 0.;
1283 hermit[hermit_dim1 + 4] = .625;
1284 hermit[hermit_dim1 + 5] = 0.;
1285 hermit[hermit_dim1 + 6] = -.1875;
1287 hermit[(hermit_dim1 << 1) + 1] = .5;
1288 hermit[(hermit_dim1 << 1) + 2] = .9375;
1289 hermit[(hermit_dim1 << 1) + 3] = 0.;
1290 hermit[(hermit_dim1 << 1) + 4] = -.625;
1291 hermit[(hermit_dim1 << 1) + 5] = 0.;
1292 hermit[(hermit_dim1 << 1) + 6] = .1875;
1294 hermit[hermit_dim1 * 3 + 1] = .3125;
1295 hermit[hermit_dim1 * 3 + 2] = -.4375;
1296 hermit[hermit_dim1 * 3 + 3] = -.375;
1297 hermit[hermit_dim1 * 3 + 4] = .625;
1298 hermit[hermit_dim1 * 3 + 5] = .0625;
1299 hermit[hermit_dim1 * 3 + 6] = -.1875;
1301 hermit[(hermit_dim1 << 2) + 1] = -.3125;
1302 hermit[(hermit_dim1 << 2) + 2] = -.4375;
1303 hermit[(hermit_dim1 << 2) + 3] = .375;
1304 hermit[(hermit_dim1 << 2) + 4] = .625;
1305 hermit[(hermit_dim1 << 2) + 5] = -.0625;
1306 hermit[(hermit_dim1 << 2) + 6] = -.1875;
1308 hermit[hermit_dim1 * 5 + 1] = .0625;
1309 hermit[hermit_dim1 * 5 + 2] = -.0625;
1310 hermit[hermit_dim1 * 5 + 3] = -.125;
1311 hermit[hermit_dim1 * 5 + 4] = .125;
1312 hermit[hermit_dim1 * 5 + 5] = .0625;
1313 hermit[hermit_dim1 * 5 + 6] = -.0625;
1315 hermit[hermit_dim1 * 6 + 1] = .0625;
1316 hermit[hermit_dim1 * 6 + 2] = .0625;
1317 hermit[hermit_dim1 * 6 + 3] = -.125;
1318 hermit[hermit_dim1 * 6 + 4] = -.125;
1319 hermit[hermit_dim1 * 6 + 5] = .0625;
1320 hermit[hermit_dim1 * 6 + 6] = .0625;
1325 /* ------------------------------ The End -------------------------------
1328 AdvApp2Var_SysBase::maermsg_("MMA1HER", iercod, 7L);
1330 AdvApp2Var_SysBase::mgsomsg_("MMA1HER", 7L);
1334 //=======================================================================
1335 //function : mma1jak_
1337 //=======================================================================
1338 int mma1jak_(integer *ndimen,
1348 /* System generated locals */
1349 integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
1350 crvjac_dim1, crvjac_offset, cgauss_dim1;
1352 /* Local variables */
1355 /* **********************************************************************
1360 /* Calculate the curve of approximation of a non-polynomial function */
1361 /* in the base of Jacobi. */
1365 /* FUNCTION,DISCRETISATION,APPROXIMATION,CONSTRAINT,CURVE,JACOBI */
1367 /* INPUT ARGUMENTS : */
1368 /* ------------------ */
1369 /* NDIMEN: Total dimension of the space (sum of dimensions */
1370 /* of sub-spaces) */
1371 /* NBROOT: Nb of points of discretization of the iso, extremities not
1373 /* IORDRE: Order of constraint at the extremities of the boundary */
1374 /* -1 = no constraints, */
1375 /* 0 = constraints of passage of limits (i.e. C0), */
1376 /* 1 = C0 + constraints of 1st derivatives (i.e. C1), */
1377 /* 2 = C1 + constraints of 2nd derivatives (i.e. C2). */
1378 /* NDGJAC: Degree of development in series to be used for calculation in the
1379 /* base of Jacobi. */
1381 /* OUTPUT ARGUMENTS : */
1382 /* ------------------- */
1383 /* CRVJAC : Curve of approximation of FONCNP with (eventually) */
1384 /* taking into account of constraints at the extremities. */
1385 /* This curve is of degree NDGJAC. */
1386 /* IERCOD : Error code : */
1387 /* 0 = All is ok. */
1388 /* 33 = Pb to return data of du block data */
1389 /* of coeff. of integration by GAUSS method. */
1390 /* by program MMAPPTT. */
1392 /* COMMONS USED : */
1393 /* ---------------- */
1395 /* REFERENCES CALLED : */
1396 /* ----------------------- */
1398 /* **********************************************************************
1401 /* Name of the routine */
1403 /* Parameter adjustments */
1404 diftab_dim1 = *nbroot / 2 + 1;
1405 diftab_offset = diftab_dim1;
1406 diftab -= diftab_offset;
1407 somtab_dim1 = *nbroot / 2 + 1;
1408 somtab_offset = somtab_dim1;
1409 somtab -= somtab_offset;
1410 crvjac_dim1 = *ndgjac + 1;
1411 crvjac_offset = crvjac_dim1;
1412 crvjac -= crvjac_offset;
1413 cgauss_dim1 = *nbroot / 2 + 1;
1416 ibb = AdvApp2Var_SysBase::mnfndeb_();
1418 AdvApp2Var_SysBase::mgenmsg_("MMA1JAK", 7L);
1422 /* ----------------- Recover coeffs of integration by Gauss -----------
1425 AdvApp2Var_ApproxF2var::mmapptt_(ndgjac, nbroot, iordre, cgauss, iercod);
1431 /* --------------- Calculate the curve in the base of Jacobi -----------
1434 mmmapcoe_(ndimen, ndgjac, iordre, nbroot, &somtab[somtab_offset], &diftab[
1435 diftab_offset], cgauss, &crvjac[crvjac_offset]);
1437 /* ------------------------------ The End -------------------------------
1442 AdvApp2Var_SysBase::maermsg_("MMA1JAK", iercod, 7L);
1445 AdvApp2Var_SysBase::mgsomsg_("MMA1JAK", 7L);
1450 //=======================================================================
1451 //function : mma1noc_
1453 //=======================================================================
1454 int mma1noc_(doublereal *dfuvin,
1463 /* System generated locals */
1468 /* Local variables */
1469 static doublereal rider, riord;
1470 static integer nd, ibb;
1471 static doublereal bid;
1472 /* **********************************************************************
1477 /* Normalization of constraints of derivatives, defined on DFUVIN */
1478 /* on block DUVOUT. */
1482 /* ALL, AB_SPECIFI::VECTEUR&,DERIVEE&,NORMALISATION,&VECTEUR */
1484 /* INPUT ARGUMENTS : */
1485 /* ------------------ */
1486 /* DFUVIN: Limits of the block of definition by U and by V where
1488 /* constraints CNTRIN are defined. */
1489 /* NDIMEN: Dimension of the space. */
1490 /* IORDRE: Order of constraint imposed at the extremities of the iso. */
1491 /* (if Iso-U, it is necessary to calculate derivatives by V and vice */
1493 /* = 0, the extremities of the iso are calculated */
1494 /* = 1, additionally the 1st derivative in the direction */
1495 /* of the iso is calculated */
1496 /* = 2, additionally the 2nd derivative in the direction */
1497 /* of the iso is calculated */
1498 /* CNTRIN: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1499 /* of order IORDRE of F(Uc,v) or of F(u,Vc), following the */
1500 /* value of ISOFAV, RENORMALIZED by u and v in (-1,1). */
1501 /* DUVOUT: Limits of the block of definition by U and by V where the */
1502 /* constraints CNTOUT will be defined. */
1503 /* ISOFAV: Isoparameter fixed for the discretization; */
1504 /* = 1, discretization with fixed U=Uc and variable V. */
1505 /* = 2, discretization with fixed V=Vc and variable U. */
1506 /* IDERIV: Ordre de derivee transverse a l'iso fixee (Si Iso-U=Uc */
1507 /* is fixed, the derivative of order IDERIV is discretized by U */
1508 /* of F(Uc,v). The same if iso-V is fixed). */
1509 /* Varies from (positioning) to 2 (2nd derivative). */
1511 /* OUTPUT ARGUMENTS : */
1512 /* ------------------- */
1513 /* CNTOUT: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1514 /* of order IORDRE of F(Uc,v) or of F(u,Vc), depending on the */
1515 /* value of ISOFAV, RENORMALIZED for u and v in DUVOUT. */
1517 /* COMMONS USED : */
1518 /* ---------------- */
1520 /* REFERENCES CALLED : */
1521 /* --------------------- */
1523 /* DESCRIPTION/NOTES/LIMITATIONS : */
1524 /* ------------------------------- */
1525 /* CNTRIN can be an output/input argument, */
1528 /* CALL MMA1NOC(DFUVIN,NDIMEN,IORDRE,CNTRIN,DUVOUT */
1529 /* 1 ,ISOFAV,IDERIV,CNTRIN) */
1533 /* **********************************************************************
1536 /* Name of the routine */
1539 /* Parameter adjustments */
1546 ibb = AdvApp2Var_SysBase::mnfndeb_();
1548 AdvApp2Var_SysBase::mgenmsg_("MMA1NOC", 7L);
1551 /* --------------- Determination of coefficients of normalization -------
1555 d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1556 rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1557 d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1558 riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1561 d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1562 rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1563 d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1564 riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1567 /* ------------- Renormalization of the vector of constraint ---------------
1570 bid = rider * riord;
1572 for (nd = 1; nd <= i__1; ++nd) {
1573 cntout[nd] = bid * cntrin[nd];
1577 /* ------------------------------ The end -------------------------------
1581 AdvApp2Var_SysBase::mgsomsg_("MMA1NOC", 7L);
1586 //=======================================================================
1587 //function : mma1nop_
1589 //=======================================================================
1590 int mma1nop_(integer *nbroot,
1598 /* System generated locals */
1601 /* Local variables */
1602 static doublereal alinu, blinu, alinv, blinv;
1603 static integer ii, ibb;
1607 /* ***********************************************************************
1612 /* Normalization of parameters of an iso, starting from */
1613 /* parametric block and parameters on (-1,1). */
1617 /* TOUS,AB_SPECIFI :: ISO&,POINT&,NORMALISATION,&POINT,&ISO */
1619 /* INPUT ARGUMENTS : */
1620 /* -------------------- */
1621 /* NBROOT: Nb of points of discretisation INSIDE the iso */
1622 /* defined on (-1,1). */
1623 /* ROOTLG: Table of discretization parameters on )-1,1( */
1625 /* UVFONC: Block of definition of the iso */
1626 /* ISOFAV: = 1, this is iso-u; =2, this is iso-v. */
1628 /* OUTPUT ARGUMENTS : */
1629 /* --------------------- */
1630 /* TTABLE: Table of parameters renormalized on UVFONC of the iso.
1632 /* IERCOD: = 0, OK */
1633 /* = 1, ISOFAV is out of allowed values. */
1636 /* **********************************************************************
1638 /* Name of the routine */
1641 /* Parameter adjustments */
1646 ibb = AdvApp2Var_SysBase::mnfndeb_();
1648 AdvApp2Var_SysBase::mgenmsg_("MMA1NOP", 7L);
1651 alinu = (uvfonc[4] - uvfonc[3]) / 2.;
1652 blinu = (uvfonc[4] + uvfonc[3]) / 2.;
1653 alinv = (uvfonc[6] - uvfonc[5]) / 2.;
1654 blinv = (uvfonc[6] + uvfonc[5]) / 2.;
1657 ttable[0] = uvfonc[5];
1659 for (ii = 1; ii <= i__1; ++ii) {
1660 ttable[ii] = alinv * rootlg[ii] + blinv;
1663 ttable[*nbroot + 1] = uvfonc[6];
1664 } else if (*isofav == 2) {
1665 ttable[0] = uvfonc[3];
1667 for (ii = 1; ii <= i__1; ++ii) {
1668 ttable[ii] = alinu * rootlg[ii] + blinu;
1671 ttable[*nbroot + 1] = uvfonc[4];
1678 /* ------------------------------ THE END -------------------------------
1687 AdvApp2Var_SysBase::maermsg_("MMA1NOP", iercod, 7L);
1690 AdvApp2Var_SysBase::mgsomsg_("MMA1NOP", 7L);
1697 //=======================================================================
1698 //function : mma2ac1_
1700 //=======================================================================
1701 int AdvApp2Var_ApproxF2var::mma2ac1_(integer const *ndimen,
1702 integer const *mxujac,
1703 integer const *mxvjac,
1704 integer const *iordru,
1705 integer const *iordrv,
1706 doublereal const *contr1,
1707 doublereal const * contr2,
1708 doublereal const *contr3,
1709 doublereal const *contr4,
1710 doublereal const *uhermt,
1711 doublereal const *vhermt,
1715 /* System generated locals */
1716 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
1717 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
1718 contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1,
1719 uhermt_offset, vhermt_dim1, vhermt_offset, patjac_dim1,
1720 patjac_dim2, patjac_offset, i__1, i__2, i__3, i__4, i__5;
1722 /* Local variables */
1723 static logical ldbg;
1724 static integer ndgu, ndgv;
1725 static doublereal bidu1, bidu2, bidv1, bidv2;
1726 static integer ioru1, iorv1, ii, nd, jj, ku, kv;
1727 static doublereal cnt1, cnt2, cnt3, cnt4;
1731 /* **********************************************************************
1736 /* Add polynoms of edge constraints. */
1740 /* TOUS,AB_SPECIFI::POINT&,CONTRAINTE&,ADDITION,&POLYNOME */
1742 /* INPUT ARGUMENTS : */
1743 /* ------------------ */
1744 /* NDIMEN: Dimension of the space. */
1745 /* MXUJAC: Max degree of the polynom of approximation by U. The */
1746 /* representation in the orthogonal base starts from degree */
1747 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1748 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1749 /* MXVJAC: Max degree of the polynom of approximation by V. The */
1750 /* representation in the orthogonal base starts from degree */
1751 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1752 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1753 /* IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
1754 /* to the step of constraints: C0, C1 or C2. */
1755 /* IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1756 /* to the step of constraints: C0, C1 or C2. */
1757 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
1758 /* extremities of F(U0,V0) and its derivatives. */
1759 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
1760 /* extremities of F(U1,V0) and its derivatives. */
1761 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
1762 /* extremities of F(U0,V1) and its derivatives. */
1763 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
1764 /* extremities of F(U1,V1) and its derivatives. */
1765 /* UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
1766 /* VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1767 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1768 /* of F(u,v) WITHOUT taking into account the constraints. */
1770 /* OUTPUT ARGUMENTS : */
1771 /* ------------------- */
1772 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1773 /* of F(u,v) WITH taking into account of constraints. */
1775 /* **********************************************************************
1777 /* Name of the routine */
1779 /* --------------------------- Initialization --------------------------
1782 /* Parameter adjustments */
1783 patjac_dim1 = *mxujac + 1;
1784 patjac_dim2 = *mxvjac + 1;
1785 patjac_offset = patjac_dim1 * patjac_dim2;
1786 patjac -= patjac_offset;
1787 uhermt_dim1 = (*iordru << 1) + 2;
1788 uhermt_offset = uhermt_dim1;
1789 uhermt -= uhermt_offset;
1790 vhermt_dim1 = (*iordrv << 1) + 2;
1791 vhermt_offset = vhermt_dim1;
1792 vhermt -= vhermt_offset;
1793 contr4_dim1 = *ndimen;
1794 contr4_dim2 = *iordru + 2;
1795 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
1796 contr4 -= contr4_offset;
1797 contr3_dim1 = *ndimen;
1798 contr3_dim2 = *iordru + 2;
1799 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
1800 contr3 -= contr3_offset;
1801 contr2_dim1 = *ndimen;
1802 contr2_dim2 = *iordru + 2;
1803 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
1804 contr2 -= contr2_offset;
1805 contr1_dim1 = *ndimen;
1806 contr1_dim2 = *iordru + 2;
1807 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
1808 contr1 -= contr1_offset;
1811 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1813 AdvApp2Var_SysBase::mgenmsg_("MMA2AC1", 7L);
1816 /* ------------ SUBTRACTION OF ANGULAR CONSTRAINTS -------------------
1819 ioru1 = *iordru + 1;
1820 iorv1 = *iordrv + 1;
1821 ndgu = (*iordru << 1) + 1;
1822 ndgv = (*iordrv << 1) + 1;
1825 for (jj = 1; jj <= i__1; ++jj) {
1827 for (ii = 1; ii <= i__2; ++ii) {
1829 for (nd = 1; nd <= i__3; ++nd) {
1830 cnt1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
1831 cnt2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
1832 cnt3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
1833 cnt4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
1835 for (kv = 0; kv <= i__4; ++kv) {
1836 bidv1 = vhermt[kv + ((jj << 1) - 1) * vhermt_dim1];
1837 bidv2 = vhermt[kv + (jj << 1) * vhermt_dim1];
1839 for (ku = 0; ku <= i__5; ++ku) {
1840 bidu1 = uhermt[ku + ((ii << 1) - 1) * uhermt_dim1];
1841 bidu2 = uhermt[ku + (ii << 1) * uhermt_dim1];
1842 patjac[ku + (kv + nd * patjac_dim2) * patjac_dim1] =
1843 patjac[ku + (kv + nd * patjac_dim2) *
1844 patjac_dim1] - bidu1 * bidv1 * cnt1 - bidu2 *
1845 bidv1 * cnt2 - bidu1 * bidv2 * cnt3 - bidu2 *
1858 /* ------------------------------ The end -------------------------------
1862 AdvApp2Var_SysBase::mgsomsg_("MMA2AC1", 7L);
1867 //=======================================================================
1868 //function : mma2ac2_
1870 //=======================================================================
1871 int AdvApp2Var_ApproxF2var::mma2ac2_(const integer *ndimen,
1872 const integer *mxujac,
1873 const integer *mxvjac,
1874 const integer *iordrv,
1875 const integer *nclimu,
1876 const integer *ncfiv1,
1877 const doublereal *crbiv1,
1878 const integer *ncfiv2,
1879 const doublereal *crbiv2,
1880 const doublereal *vhermt,
1884 /* System generated locals */
1885 integer crbiv1_dim1, crbiv1_dim2, crbiv1_offset, crbiv2_dim1, crbiv2_dim2,
1886 crbiv2_offset, patjac_dim1, patjac_dim2, patjac_offset,
1887 vhermt_dim1, vhermt_offset, i__1, i__2, i__3, i__4;
1889 /* Local variables */
1890 static logical ldbg;
1891 static integer ndgv1, ndgv2, ii, jj, nd, kk;
1892 static doublereal bid1, bid2;
1894 /* **********************************************************************
1899 /* Add polynoms of constraints */
1903 /* FUNCTION,APPROXIMATION,COEFFICIENT,POLYNOM */
1905 /* INPUT ARGUMENTS : */
1906 /* ------------------ */
1907 /* NDIMEN: Dimension of the space. */
1908 /* MXUJAC: Max degree of the polynom of approximation by U. The */
1909 /* representation in the orthogonal base starts from degree */
1910 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1911 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1912 /* MXVJAC: Max degree of the polynom of approximation by V. The */
1913 /* representation in the orthogonal base starts from degree */
1914 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1915 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1916 /* IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1917 /* to the step of constraints: C0, C1 or C2. */
1918 /* NCLIMU LIMIT nb of coeff by u of the solution P(u,v)
1919 * NCFIV1: Nb of Coeff. of curves stored in CRBIV1. */
1920 /* CRBIV1: Table of coeffs of the approximation of iso-V0 and its */
1921 /* derivatives till order IORDRV. */
1922 /* NCFIV2: Nb of Coeff. of curves stored in CRBIV2. */
1923 /* CRBIV2: Table of coeffs of approximation of iso-V1 and its */
1924 /* derivatives till order IORDRV. */
1925 /* VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1926 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1927 /* of F(u,v) WITHOUT taking into account the constraints. */
1929 /* OUTPUT ARGUMENTS : */
1930 /* ------------------- */
1931 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1932 /* of F(u,v) WITH taking into account of constraints. */
1937 /* **********************************************************************
1939 /* Name of the routine */
1941 /* --------------------------- Initialisations --------------------------
1944 /* Parameter adjustments */
1945 patjac_dim1 = *mxujac + 1;
1946 patjac_dim2 = *mxvjac + 1;
1947 patjac_offset = patjac_dim1 * patjac_dim2;
1948 patjac -= patjac_offset;
1949 vhermt_dim1 = (*iordrv << 1) + 2;
1950 vhermt_offset = vhermt_dim1;
1951 vhermt -= vhermt_offset;
1954 crbiv2_dim1 = *nclimu;
1955 crbiv2_dim2 = *ndimen;
1956 crbiv2_offset = crbiv2_dim1 * (crbiv2_dim2 + 1);
1957 crbiv2 -= crbiv2_offset;
1958 crbiv1_dim1 = *nclimu;
1959 crbiv1_dim2 = *ndimen;
1960 crbiv1_offset = crbiv1_dim1 * (crbiv1_dim2 + 1);
1961 crbiv1 -= crbiv1_offset;
1964 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1966 AdvApp2Var_SysBase::mgenmsg_("MMA2AC2", 7L);
1969 /* ------------ ADDING of coeff by u of curves, by v of Hermit --------
1973 for (ii = 1; ii <= i__1; ++ii) {
1974 ndgv1 = ncfiv1[ii] - 1;
1975 ndgv2 = ncfiv2[ii] - 1;
1977 for (nd = 1; nd <= i__2; ++nd) {
1978 i__3 = (*iordrv << 1) + 1;
1979 for (jj = 0; jj <= i__3; ++jj) {
1980 bid1 = vhermt[jj + ((ii << 1) - 1) * vhermt_dim1];
1982 for (kk = 0; kk <= i__4; ++kk) {
1983 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
1984 bid1 * crbiv1[kk + (nd + ii * crbiv1_dim2) *
1988 bid2 = vhermt[jj + (ii << 1) * vhermt_dim1];
1990 for (kk = 0; kk <= i__4; ++kk) {
1991 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
1992 bid2 * crbiv2[kk + (nd + ii * crbiv2_dim2) *
2003 /* ------------------------------ The end -------------------------------
2007 AdvApp2Var_SysBase::mgsomsg_("MMA2AC2", 7L);
2013 //=======================================================================
2014 //function : mma2ac3_
2016 //=======================================================================
2017 int AdvApp2Var_ApproxF2var::mma2ac3_(const integer *ndimen,
2018 const integer *mxujac,
2019 const integer *mxvjac,
2020 const integer *iordru,
2021 const integer *nclimv,
2022 const integer *ncfiu1,
2023 const doublereal * crbiu1,
2024 const integer *ncfiu2,
2025 const doublereal *crbiu2,
2026 const doublereal *uhermt,
2030 /* System generated locals */
2031 integer crbiu1_dim1, crbiu1_dim2, crbiu1_offset, crbiu2_dim1, crbiu2_dim2,
2032 crbiu2_offset, patjac_dim1, patjac_dim2, patjac_offset,
2033 uhermt_dim1, uhermt_offset, i__1, i__2, i__3, i__4;
2035 /* Local variables */
2036 static logical ldbg;
2037 static integer ndgu1, ndgu2, ii, jj, nd, kk;
2038 static doublereal bid1, bid2;
2043 /* **********************************************************************
2048 /* Ajout des polynomes de contraintes */
2052 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
2054 /* INPUT ARGUMENTS : */
2055 /* ------------------ */
2056 /* NDIMEN: Dimension of the space. */
2057 /* MXUJAC: Max degree of the polynom of approximation by U. The */
2058 /* representation in the orthogonal base starts from degree */
2059 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2060 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2061 /* MXVJAC: Max degree of the polynom of approximation by V. The */
2062 /* representation in the orthogonal base starts from degree */
2063 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2064 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2065 /* IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
2066 /* to the step of constraints: C0, C1 or C2. */
2067 /* NCLIMV LIMIT nb of coeff by v of the solution P(u,v)
2068 * NCFIU1: Nb of Coeff. of curves stored in CRBIU1. */
2069 /* CRBIU1: Table of coeffs of the approximation of iso-U0 and its */
2070 /* derivatives till order IORDRU. */
2071 /* NCFIU2: Nb of Coeff. of curves stored in CRBIU2. */
2072 /* CRBIU2: Table of coeffs of approximation of iso-U1 and its */
2073 /* derivatives till order IORDRU */
2074 /* UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
2075 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
2076 /* of F(u,v) WITHOUT taking into account the constraints. */
2078 /* OUTPUT ARGUMENTS : */
2079 /* ------------------- */
2080 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
2081 /* of F(u,v) WITH taking into account of constraints. */
2085 /* **********************************************************************
2087 /* The name of the routine */
2089 /* --------------------------- Initializations --------------------------
2092 /* Parameter adjustments */
2093 patjac_dim1 = *mxujac + 1;
2094 patjac_dim2 = *mxvjac + 1;
2095 patjac_offset = patjac_dim1 * patjac_dim2;
2096 patjac -= patjac_offset;
2097 uhermt_dim1 = (*iordru << 1) + 2;
2098 uhermt_offset = uhermt_dim1;
2099 uhermt -= uhermt_offset;
2102 crbiu2_dim1 = *nclimv;
2103 crbiu2_dim2 = *ndimen;
2104 crbiu2_offset = crbiu2_dim1 * (crbiu2_dim2 + 1);
2105 crbiu2 -= crbiu2_offset;
2106 crbiu1_dim1 = *nclimv;
2107 crbiu1_dim2 = *ndimen;
2108 crbiu1_offset = crbiu1_dim1 * (crbiu1_dim2 + 1);
2109 crbiu1 -= crbiu1_offset;
2112 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
2114 AdvApp2Var_SysBase::mgenmsg_("MMA2AC3", 7L);
2117 /* ------------ ADDING of coeff by u of curves, by v of Hermit --------
2121 for (ii = 1; ii <= i__1; ++ii) {
2122 ndgu1 = ncfiu1[ii] - 1;
2123 ndgu2 = ncfiu2[ii] - 1;
2125 for (nd = 1; nd <= i__2; ++nd) {
2127 for (jj = 0; jj <= i__3; ++jj) {
2128 bid1 = crbiu1[jj + (nd + ii * crbiu1_dim2) * crbiu1_dim1];
2129 i__4 = (*iordru << 1) + 1;
2130 for (kk = 0; kk <= i__4; ++kk) {
2131 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
2132 bid1 * uhermt[kk + ((ii << 1) - 1) * uhermt_dim1];
2138 for (jj = 0; jj <= i__3; ++jj) {
2139 bid2 = crbiu2[jj + (nd + ii * crbiu2_dim2) * crbiu2_dim1];
2140 i__4 = (*iordru << 1) + 1;
2141 for (kk = 0; kk <= i__4; ++kk) {
2142 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
2143 bid2 * uhermt[kk + (ii << 1) * uhermt_dim1];
2154 /* ------------------------------ The end -------------------------------
2158 AdvApp2Var_SysBase::mgsomsg_("MMA2AC3", 7L);
2163 //=======================================================================
2164 //function : mma2can_
2166 //=======================================================================
2167 int AdvApp2Var_ApproxF2var::mma2can_(const integer *ncfmxu,
2168 const integer *ncfmxv,
2169 const integer *ndimen,
2170 const integer *iordru,
2171 const integer *iordrv,
2172 const integer *ncoefu,
2173 const integer *ncoefv,
2174 const doublereal *patjac,
2180 /* System generated locals */
2181 integer patjac_dim1, patjac_dim2, patjac_offset, patcan_dim1, patcan_dim2,
2182 patcan_offset, i__1, i__2;
2184 /* Local variables */
2185 static logical ldbg;
2186 static integer ilon1, ilon2, ii, nd;
2191 /* **********************************************************************
2196 /* Change of Jacobi base to canonical (-1,1) and writing in a greater */
2201 /* ALL,AB_SPECIFI,CARREAU&,CONVERSION,JACOBI,CANNONIQUE,&CARREAU */
2203 /* INPUT ARGUMENTS : */
2204 /* -------------------- */
2205 /* NCFMXU: Dimension by U of resulting table PATCAN */
2206 /* NCFMXV: Dimension by V of resulting table PATCAN */
2207 /* NDIMEN: Dimension of the workspace. */
2208 /* IORDRU: Order of constraint by U */
2209 /* IORDRV: Order of constraint by V. */
2210 /* NCOEFU: Nb of coeff by U of square PATJAC */
2211 /* NCOEFV: Nb of coeff by V of square PATJAC */
2212 /* PATJAC: Square in the base of Jacobi of order IORDRU by U and */
2215 /* OUTPUT ARGUMENTS : */
2216 /* --------------------- */
2217 /* PATAUX: Auxiliary Table. */
2218 /* PATCAN: Table of coefficients in the canonic base. */
2219 /* IERCOD: Error code. */
2220 /* = 0, everything goes well, and all things are equal. */
2221 /* = 1, the program refuses to process with incorrect input arguments */
2224 /* COMMONS USED : */
2225 /* ------------------ */
2227 /* REFERENCES CALLED : */
2228 /* --------------------- */
2230 /* DESCRIPTION/NOTES/LIMITATIONS : */
2231 /* ----------------------------------- */
2233 /* **********************************************************************
2237 /* Parameter adjustments */
2238 patcan_dim1 = *ncfmxu;
2239 patcan_dim2 = *ncfmxv;
2240 patcan_offset = patcan_dim1 * (patcan_dim2 + 1) + 1;
2241 patcan -= patcan_offset;
2243 patjac_dim1 = *ncoefu;
2244 patjac_dim2 = *ncoefv;
2245 patjac_offset = patjac_dim1 * (patjac_dim2 + 1) + 1;
2246 patjac -= patjac_offset;
2249 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
2251 AdvApp2Var_SysBase::mgenmsg_("MMA2CAN", 7L);
2255 if (*iordru < -1 || *iordru > 2) {
2258 if (*iordrv < -1 || *iordrv > 2) {
2261 if (*ncoefu > *ncfmxu || *ncoefv > *ncfmxv) {
2265 /* --> Pass to canonic base (-1,1) */
2266 mmjacpt_(ndimen, ncoefu, ncoefv, iordru, iordrv, &patjac[patjac_offset], &
2267 pataux[1], &patcan[patcan_offset]);
2269 /* --> Write all in a greater table */
2270 AdvApp2Var_MathBase::mmfmca8_((integer *)ncoefu,
2276 (doublereal *)&patcan[patcan_offset],
2277 (doublereal *)&patcan[patcan_offset]);
2279 /* --> Complete with zeros the resulting table. */
2280 ilon1 = *ncfmxu - *ncoefu;
2281 ilon2 = *ncfmxu * (*ncfmxv - *ncoefv);
2283 for (nd = 1; nd <= i__1; ++nd) {
2286 for (ii = 1; ii <= i__2; ++ii) {
2287 AdvApp2Var_SysBase::mvriraz_(&ilon1,
2288 (char *)&patcan[*ncoefu + 1 + (ii + nd * patcan_dim2) * patcan_dim1]);
2293 AdvApp2Var_SysBase::mvriraz_(&ilon2,
2294 (char *)&patcan[(*ncoefv + 1 + nd * patcan_dim2) * patcan_dim1 + 1]);
2301 /* ----------------------
2309 AdvApp2Var_SysBase::maermsg_("MMA2CAN", iercod, 7L);
2311 AdvApp2Var_SysBase::mgsomsg_("MMA2CAN", 7L);
2316 //=======================================================================
2317 //function : mma2cd1_
2319 //=======================================================================
2320 int mma2cd1_(integer *ndimen,
2341 static integer c__1 = 1;
2343 /* System generated locals */
2344 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
2345 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
2346 contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1,
2347 uhermt_offset, vhermt_dim1, vhermt_offset, fpntbu_dim1,
2348 fpntbu_offset, fpntbv_dim1, fpntbv_offset, sosotb_dim1,
2349 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2350 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2351 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4,
2354 /* Local variables */
2355 static integer ncfhu, ncfhv, nuroo, nvroo, nd, ii, jj, kk, ll, ibb, kkm,
2357 static doublereal bid1, bid2, bid3, bid4;
2358 static doublereal diu1, diu2, div1, div2, sou1, sou2, sov1, sov2;
2363 /* **********************************************************************
2368 /* Discretisation on the parameters of polynoms of interpolation */
2369 /* of constraints at the corners of order IORDRE. */
2373 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2375 /* INPUT ARGUMENTS : */
2376 /* ------------------ */
2377 /* NDIMEN: Dimension of the space. */
2378 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2379 /* This is also the nb of root of Legendre polynom where discretization is done. */
2380 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
2382 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
2383 /* This is also the nb of root of Legendre polynom where discretization is done. */
2384 /* VROOTL: Table of discretization parameters on (-1,1) by V.
2385 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
2386 /* = 0, calculate the extremities of iso-V */
2387 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2388 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2389 /* IORDRV: Order of constraint imposed at the extremities of iso-U */
2390 /* = 0, calculate the extremities of iso-U */
2391 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
2392 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
2393 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
2394 /* extremities of F(U0,V0) and its derivatives. */
2395 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
2396 /* extremities of F(U1,V0) and its derivatives. */
2397 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
2398 /* extremities of F(U0,V1) and its derivatives. */
2399 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
2400 /* extremities of F(U1,V1) and its derivatives. */
2401 /* SOSOTB: Preinitialized table (input/output argument). */
2402 /* DISOTB: Preinitialized table (input/output argument). */
2403 /* SODITB: Preinitialized table (input/output argument). */
2404 /* DIDITB: Preinitialized table (input/output argument) */
2406 /* OUTPUT ARGUMENTS : */
2407 /* ------------------- */
2408 /* FPNTBU: Auxiliary table. */
2409 /* FPNTBV: Auxiliary table. */
2410 /* UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
2411 /* VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2412 /* SOSOTB: Table where the terms of constraints are added */
2413 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2414 /* with ui and vj positive roots of the Legendre polynom */
2415 /* of degree NBPNTU and NBPNTV respectively. */
2416 /* DISOTB: Table where the terms of constraints are added */
2417 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2418 /* with ui and vj positive roots of the polynom of Legendre */
2419 /* of degree NBPNTU and NBPNTV respectively. */
2420 /* SODITB: Table where the terms of constraints are added */
2421 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2422 /* with ui and vj positive roots of the polynom of Legendre */
2423 /* of degree NBPNTU and NBPNTV respectively. */
2424 /* DIDITB: Table where the terms of constraints are added */
2425 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2426 /* with ui and vj positive roots of the polynom of Legendre */
2427 /* of degree NBPNTU and NBPNTV respectively. */
2429 /* COMMONS USED : */
2430 /* ---------------- */
2432 /* REFERENCES CALLED : */
2433 /* ----------------------- */
2435 /* DESCRIPTION/NOTES/LIMITATIONS : */
2436 /* ----------------------------------- */
2439 /* **********************************************************************
2442 /* Name of the routine */
2445 /* Parameter adjustments */
2447 diditb_dim1 = *nbpntu / 2 + 1;
2448 diditb_dim2 = *nbpntv / 2 + 1;
2449 diditb_offset = diditb_dim1 * diditb_dim2;
2450 diditb -= diditb_offset;
2451 disotb_dim1 = *nbpntu / 2;
2452 disotb_dim2 = *nbpntv / 2;
2453 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2454 disotb -= disotb_offset;
2455 soditb_dim1 = *nbpntu / 2;
2456 soditb_dim2 = *nbpntv / 2;
2457 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2458 soditb -= soditb_offset;
2459 sosotb_dim1 = *nbpntu / 2 + 1;
2460 sosotb_dim2 = *nbpntv / 2 + 1;
2461 sosotb_offset = sosotb_dim1 * sosotb_dim2;
2462 sosotb -= sosotb_offset;
2464 uhermt_dim1 = (*iordru << 1) + 2;
2465 uhermt_offset = uhermt_dim1;
2466 uhermt -= uhermt_offset;
2467 fpntbu_dim1 = *nbpntu;
2468 fpntbu_offset = fpntbu_dim1 + 1;
2469 fpntbu -= fpntbu_offset;
2470 vhermt_dim1 = (*iordrv << 1) + 2;
2471 vhermt_offset = vhermt_dim1;
2472 vhermt -= vhermt_offset;
2473 fpntbv_dim1 = *nbpntv;
2474 fpntbv_offset = fpntbv_dim1 + 1;
2475 fpntbv -= fpntbv_offset;
2476 contr4_dim1 = *ndimen;
2477 contr4_dim2 = *iordru + 2;
2478 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
2479 contr4 -= contr4_offset;
2480 contr3_dim1 = *ndimen;
2481 contr3_dim2 = *iordru + 2;
2482 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
2483 contr3 -= contr3_offset;
2484 contr2_dim1 = *ndimen;
2485 contr2_dim2 = *iordru + 2;
2486 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
2487 contr2 -= contr2_offset;
2488 contr1_dim1 = *ndimen;
2489 contr1_dim2 = *iordru + 2;
2490 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
2491 contr1 -= contr1_offset;
2494 ibb = AdvApp2Var_SysBase::mnfndeb_();
2496 AdvApp2Var_SysBase::mgenmsg_("MMA2CD1", 7L);
2499 /* ------------------- Discretisation of Hermite polynoms -----------
2502 ncfhu = (*iordru + 1) << 1;
2504 for (ii = 1; ii <= i__1; ++ii) {
2506 for (ll = 1; ll <= i__2; ++ll) {
2507 AdvApp2Var_MathBase::mmmpocur_(&ncfhu, &c__1, &ncfhu, &uhermt[ii * uhermt_dim1], &
2508 urootl[ll], &fpntbu[ll + ii * fpntbu_dim1]);
2513 ncfhv = (*iordrv + 1) << 1;
2515 for (jj = 1; jj <= i__1; ++jj) {
2517 for (kk = 1; kk <= i__2; ++kk) {
2518 AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[jj * vhermt_dim1], &
2519 vrootl[kk], &fpntbv[kk + jj * fpntbv_dim1]);
2525 /* ---- The discretizations of polynoms of constraints are subtracted ----
2528 nuroo = *nbpntu / 2;
2529 nvroo = *nbpntv / 2;
2531 for (nd = 1; nd <= i__1; ++nd) {
2534 for (jj = 1; jj <= i__2; ++jj) {
2536 for (ii = 1; ii <= i__3; ++ii) {
2537 bid1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
2538 bid2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
2539 bid3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
2540 bid4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
2543 for (kk = 1; kk <= i__4; ++kk) {
2544 kkp = (*nbpntv + 1) / 2 + kk;
2545 kkm = nvroo - kk + 1;
2546 sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] +
2547 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2548 div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] -
2549 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2550 sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[kkm
2551 + (jj << 1) * fpntbv_dim1];
2552 div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[kkm
2553 + (jj << 1) * fpntbv_dim1];
2555 for (ll = 1; ll <= i__5; ++ll) {
2556 llp = (*nbpntu + 1) / 2 + ll;
2557 llm = nuroo - ll + 1;
2558 sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] +
2559 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2560 diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] -
2561 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2562 sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2563 llm + (ii << 1) * fpntbu_dim1];
2564 diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2565 llm + (ii << 1) * fpntbu_dim1];
2566 sosotb[ll + (kk + nd * sosotb_dim2) * sosotb_dim1] =
2567 sosotb[ll + (kk + nd * sosotb_dim2) *
2568 sosotb_dim1] - bid1 * sou1 * sov1 - bid2 *
2569 sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2571 soditb[ll + (kk + nd * soditb_dim2) * soditb_dim1] =
2572 soditb[ll + (kk + nd * soditb_dim2) *
2573 soditb_dim1] - bid1 * sou1 * div1 - bid2 *
2574 sou2 * div1 - bid3 * sou1 * div2 - bid4 *
2576 disotb[ll + (kk + nd * disotb_dim2) * disotb_dim1] =
2577 disotb[ll + (kk + nd * disotb_dim2) *
2578 disotb_dim1] - bid1 * diu1 * sov1 - bid2 *
2579 diu2 * sov1 - bid3 * diu1 * sov2 - bid4 *
2581 diditb[ll + (kk + nd * diditb_dim2) * diditb_dim1] =
2582 diditb[ll + (kk + nd * diditb_dim2) *
2583 diditb_dim1] - bid1 * diu1 * div1 - bid2 *
2584 diu2 * div1 - bid3 * diu1 * div2 - bid4 *
2591 /* ------------ Case when the discretization is done only on the roots
2593 /* ---------- of Legendre polynom of uneven degree, 0 is root
2596 if (*nbpntu % 2 == 1) {
2597 sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2598 sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2600 for (kk = 1; kk <= i__4; ++kk) {
2601 kkp = (*nbpntv + 1) / 2 + kk;
2602 kkm = nvroo - kk + 1;
2603 sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] +
2604 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2605 div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] -
2606 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2607 sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[
2608 kkm + (jj << 1) * fpntbv_dim1];
2609 div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[
2610 kkm + (jj << 1) * fpntbv_dim1];
2611 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1] =
2612 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1]
2613 - bid1 * sou1 * sov1 - bid2 * sou2 * sov1 -
2614 bid3 * sou1 * sov2 - bid4 * sou2 * sov2;
2615 diditb[(kk + nd * diditb_dim2) * diditb_dim1] =
2616 diditb[(kk + nd * diditb_dim2) * diditb_dim1]
2617 - bid1 * sou1 * div1 - bid2 * sou2 * div1 -
2618 bid3 * sou1 * div2 - bid4 * sou2 * div2;
2623 if (*nbpntv % 2 == 1) {
2624 sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2625 sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2627 for (ll = 1; ll <= i__4; ++ll) {
2628 llp = (*nbpntu + 1) / 2 + ll;
2629 llm = nuroo - ll + 1;
2630 sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] +
2631 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2632 diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] -
2633 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2634 sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2635 llm + (ii << 1) * fpntbu_dim1];
2636 diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2637 llm + (ii << 1) * fpntbu_dim1];
2638 sosotb[ll + nd * sosotb_dim2 * sosotb_dim1] = sosotb[
2639 ll + nd * sosotb_dim2 * sosotb_dim1] - bid1 *
2640 sou1 * sov1 - bid2 * sou2 * sov1 - bid3 *
2641 sou1 * sov2 - bid4 * sou2 * sov2;
2642 diditb[ll + nd * diditb_dim2 * diditb_dim1] = diditb[
2643 ll + nd * diditb_dim2 * diditb_dim1] - bid1 *
2644 diu1 * sov1 - bid2 * diu2 * sov1 - bid3 *
2645 diu1 * sov2 - bid4 * diu2 * sov2;
2650 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2651 sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2652 sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2653 sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2654 sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2655 sosotb[nd * sosotb_dim2 * sosotb_dim1] = sosotb[nd *
2656 sosotb_dim2 * sosotb_dim1] - bid1 * sou1 * sov1 -
2657 bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2659 diditb[nd * diditb_dim2 * diditb_dim1] = diditb[nd *
2660 diditb_dim2 * diditb_dim1] - bid1 * sou1 * sov1 -
2661 bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2673 /* ------------------------------ The End -------------------------------
2678 AdvApp2Var_SysBase::mgsomsg_("MMA2CD1", 7L);
2683 //=======================================================================
2684 //function : mma2cd2_
2686 //=======================================================================
2687 int mma2cd2_(integer *ndimen,
2704 static integer c__1 = 1;
2705 /* System generated locals */
2706 integer sotbv1_dim1, sotbv1_dim2, sotbv1_offset, sotbv2_dim1, sotbv2_dim2,
2707 sotbv2_offset, ditbv1_dim1, ditbv1_dim2, ditbv1_offset,
2708 ditbv2_dim1, ditbv2_dim2, ditbv2_offset, fpntab_dim1,
2709 fpntab_offset, vhermt_dim1, vhermt_offset, sosotb_dim1,
2710 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2711 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2712 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2714 /* Local variables */
2715 static integer ncfhv, nuroo, nvroo, ii, nd, jj, kk, ibb, jjm, jjp;
2716 static doublereal bid1, bid2, bid3, bid4;
2718 /* **********************************************************************
2722 /* Discretisation on the parameters of polynoms of interpolation */
2723 /* of constraints on 2 borders iso-V of order IORDRV. */
2728 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2732 /* INPUT ARGUMENTS : */
2733 /* ------------------ */
2734 /* NDIMEN: Dimension of the space. */
2735 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2736 /* This is also the nb of root of Legendre polynom where discretization is done. */
2737 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
2739 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
2740 /* This is also the nb of root of Legendre polynom where discretization is done. */
2741 /* VROOTL: Table of discretization parameters on (-1,1) by V.
2742 /* IORDRV: Order of constraint imposed at the extremities of iso-V */
2743 /* = 0, calculate the extremities of iso-V */
2744 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2745 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2746 /* SOTBV1: Table of NBPNTV/2 sums of 2 index points */
2747 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2748 /* SOTBV2: Table of NBPNTV/2 sums of 2 index points */
2749 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2750 /* DITBV1: Table of NBPNTV/2 differences of 2 index points */
2751 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2752 /* DITBV2: Table of NBPNTV/2 differences of 2 index points */
2753 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2754 /* SOSOTB: Preinitialized table (input/output argument). */
2755 /* DISOTB: Preinitialized table (input/output argument). */
2756 /* SODITB: Preinitialized table (input/output argument). */
2757 /* DIDITB: Preinitialized table (input/output argument) */
2759 /* OUTPUT ARGUMENTS : */
2760 /* ------------------- */
2761 /* FPNTAB: Auxiliary table. */
2762 /* VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2763 /* SOSOTB: Table where the terms of constraints are added */
2764 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2765 /* with ui and vj positive roots of the Legendre polynom */
2766 /* of degree NBPNTU and NBPNTV respectively. */
2767 /* DISOTB: Table where the terms of constraints are added */
2768 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2769 /* with ui and vj positive roots of the polynom of Legendre */
2770 /* of degree NBPNTU and NBPNTV respectively. */
2771 /* SODITB: Table where the terms of constraints are added */
2772 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2773 /* with ui and vj positive roots of the polynom of Legendre */
2774 /* of degree NBPNTU and NBPNTV respectively. */
2775 /* DIDITB: Table where the terms of constraints are added */
2776 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2777 /* with ui and vj positive roots of the polynom of Legendre */
2778 /* of degree NBPNTU and NBPNTV respectively. */
2780 /* COMMONS USED : */
2781 /* ---------------- */
2783 /* REFERENCES CALLED : */
2784 /* ----------------------- */
2786 /* DESCRIPTION/NOTES/LIMITATIONS : */
2787 /* ----------------------------------- */
2791 /* **********************************************************************
2794 /* Name of the routine */
2797 /* Parameter adjustments */
2798 diditb_dim1 = *nbpntu / 2 + 1;
2799 diditb_dim2 = *nbpntv / 2 + 1;
2800 diditb_offset = diditb_dim1 * diditb_dim2;
2801 diditb -= diditb_offset;
2802 disotb_dim1 = *nbpntu / 2;
2803 disotb_dim2 = *nbpntv / 2;
2804 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2805 disotb -= disotb_offset;
2806 soditb_dim1 = *nbpntu / 2;
2807 soditb_dim2 = *nbpntv / 2;
2808 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2809 soditb -= soditb_offset;
2810 sosotb_dim1 = *nbpntu / 2 + 1;
2811 sosotb_dim2 = *nbpntv / 2 + 1;
2812 sosotb_offset = sosotb_dim1 * sosotb_dim2;
2813 sosotb -= sosotb_offset;
2815 vhermt_dim1 = (*iordrv << 1) + 2;
2816 vhermt_offset = vhermt_dim1;
2817 vhermt -= vhermt_offset;
2818 fpntab_dim1 = *nbpntv;
2819 fpntab_offset = fpntab_dim1 + 1;
2820 fpntab -= fpntab_offset;
2821 ditbv2_dim1 = *nbpntu / 2 + 1;
2822 ditbv2_dim2 = *ndimen;
2823 ditbv2_offset = ditbv2_dim1 * (ditbv2_dim2 + 1);
2824 ditbv2 -= ditbv2_offset;
2825 ditbv1_dim1 = *nbpntu / 2 + 1;
2826 ditbv1_dim2 = *ndimen;
2827 ditbv1_offset = ditbv1_dim1 * (ditbv1_dim2 + 1);
2828 ditbv1 -= ditbv1_offset;
2829 sotbv2_dim1 = *nbpntu / 2 + 1;
2830 sotbv2_dim2 = *ndimen;
2831 sotbv2_offset = sotbv2_dim1 * (sotbv2_dim2 + 1);
2832 sotbv2 -= sotbv2_offset;
2833 sotbv1_dim1 = *nbpntu / 2 + 1;
2834 sotbv1_dim2 = *ndimen;
2835 sotbv1_offset = sotbv1_dim1 * (sotbv1_dim2 + 1);
2836 sotbv1 -= sotbv1_offset;
2839 ibb = AdvApp2Var_SysBase::mnfndeb_();
2841 AdvApp2Var_SysBase::mgenmsg_("MMA2CD2", 7L);
2844 /* ------------------- Discretization of Hermit polynoms -----------
2847 ncfhv = (*iordrv + 1) << 1;
2849 for (ii = 1; ii <= i__1; ++ii) {
2851 for (jj = 1; jj <= i__2; ++jj) {
2852 AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[ii * vhermt_dim1], &
2853 vrootl[jj], &fpntab[jj + ii * fpntab_dim1]);
2859 /* ---- The discretizations of polynoms of constraints are subtracted ----
2862 nuroo = *nbpntu / 2;
2863 nvroo = *nbpntv / 2;
2866 for (nd = 1; nd <= i__1; ++nd) {
2868 for (ii = 1; ii <= i__2; ++ii) {
2871 for (kk = 1; kk <= i__3; ++kk) {
2872 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1];
2873 bid2 = sotbv2[kk + (nd + ii * sotbv2_dim2) * sotbv2_dim1];
2874 bid3 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1];
2875 bid4 = ditbv2[kk + (nd + ii * ditbv2_dim2) * ditbv2_dim1];
2877 for (jj = 1; jj <= i__4; ++jj) {
2878 jjp = (*nbpntv + 1) / 2 + jj;
2879 jjm = nvroo - jj + 1;
2880 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] =
2881 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
2882 - bid1 * (fpntab[jjp + ((ii << 1) - 1) *
2883 fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) *
2884 fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) *
2885 fpntab_dim1] + fpntab[jjm + (ii << 1) *
2887 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] =
2888 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
2889 - bid3 * (fpntab[jjp + ((ii << 1) - 1) *
2890 fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) *
2891 fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) *
2892 fpntab_dim1] + fpntab[jjm + (ii << 1) *
2894 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] =
2895 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
2896 - bid1 * (fpntab[jjp + ((ii << 1) - 1) *
2897 fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) *
2898 fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) *
2899 fpntab_dim1] - fpntab[jjm + (ii << 1) *
2901 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] =
2902 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
2903 - bid3 * (fpntab[jjp + ((ii << 1) - 1) *
2904 fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) *
2905 fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) *
2906 fpntab_dim1] - fpntab[jjm + (ii << 1) *
2915 /* ------------ Case when the discretization is done only on the roots */
2916 /* ---------- of Legendre polynom of uneven degree, 0 is root */
2919 if (*nbpntv % 2 == 1) {
2921 for (ii = 1; ii <= i__2; ++ii) {
2923 for (kk = 1; kk <= i__3; ++kk) {
2924 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1]
2925 * fpntab[nvroo + 1 + ((ii << 1) - 1) *
2926 fpntab_dim1] + sotbv2[kk + (nd + ii * sotbv2_dim2)
2927 * sotbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) *
2929 sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2930 bid2 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1]
2931 * fpntab[nvroo + 1 + ((ii << 1) - 1) *
2932 fpntab_dim1] + ditbv2[kk + (nd + ii * ditbv2_dim2)
2933 * ditbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) *
2935 diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
2942 if (*nbpntu % 2 == 1) {
2944 for (ii = 1; ii <= i__2; ++ii) {
2946 for (jj = 1; jj <= i__3; ++jj) {
2947 jjp = (*nbpntv + 1) / 2 + jj;
2948 jjm = nvroo - jj + 1;
2949 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2950 fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] +
2951 fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) +
2952 sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2953 fpntab[jjp + (ii << 1) * fpntab_dim1] + fpntab[
2954 jjm + (ii << 1) * fpntab_dim1]);
2955 sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
2956 bid2 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2957 fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] -
2958 fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) +
2959 sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2960 fpntab[jjp + (ii << 1) * fpntab_dim1] - fpntab[
2961 jjm + (ii << 1) * fpntab_dim1]);
2962 diditb[jj + nd * diditb_dim2 * diditb_dim1] -= bid2;
2969 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2971 for (ii = 1; ii <= i__2; ++ii) {
2972 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * fpntab[
2973 nvroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbv2[(
2974 nd + ii * sotbv2_dim2) * sotbv2_dim1] * fpntab[nvroo
2975 + 1 + (ii << 1) * fpntab_dim1];
2976 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2985 /* ------------------------------ The End -------------------------------
2990 AdvApp2Var_SysBase::mgsomsg_("MMA2CD2", 7L);
2995 //=======================================================================
2996 //function : mma2cd3_
2998 //=======================================================================
2999 int mma2cd3_(integer *ndimen,
3016 static integer c__1 = 1;
3018 /* System generated locals */
3019 integer sotbu1_dim1, sotbu1_dim2, sotbu1_offset, sotbu2_dim1, sotbu2_dim2,
3020 sotbu2_offset, ditbu1_dim1, ditbu1_dim2, ditbu1_offset,
3021 ditbu2_dim1, ditbu2_dim2, ditbu2_offset, fpntab_dim1,
3022 fpntab_offset, uhermt_dim1, uhermt_offset, sosotb_dim1,
3023 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
3024 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3025 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
3027 /* Local variables */
3028 static integer ncfhu, nuroo, nvroo, ii, nd, jj, kk, ibb, kkm, kkp;
3029 static doublereal bid1, bid2, bid3, bid4;
3031 /* **********************************************************************
3035 /* Discretisation on the parameters of polynoms of interpolation */
3036 /* of constraints on 2 borders iso-U of order IORDRU. */
3041 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3043 /* INPUT ARGUMENTS : */
3044 /* ------------------ */
3045 /* NDIMEN: Dimension of the space. */
3046 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3047 /* This is also the nb of root of Legendre polynom where discretization is done. */
3048 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
3050 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
3051 /* This is also the nb of root of Legendre polynom where discretization is done. */
3052 /* IORDRV: Order of constraint imposed at the extremities of iso-V */
3053 /* = 0, calculate the extremities of iso-V */
3054 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3055 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3056 /* SOTBU1: Table of NBPNTU/2 sums of 2 index points */
3057 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3058 /* SOTBU2: Table of NBPNTV/2 sums of 2 index points */
3059 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3060 /* DITBU1: Table of NBPNTU/2 differences of 2 index points */
3061 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3062 /* DITBU2: Table of NBPNTU/2 differences of 2 index points */
3063 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3064 /* SOSOTB: Preinitialized table (input/output argument). */
3065 /* DISOTB: Preinitialized table (input/output argument). */
3066 /* SODITB: Preinitialized table (input/output argument). */
3067 /* DIDITB: Preinitialized table (input/output argument) */
3069 /* OUTPUT ARGUMENTS : */
3070 /* ------------------- */
3071 /* FPNTAB: Auxiliary table. */
3072 /* UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
3073 /* SOSOTB: Table where the terms of constraints are added */
3074 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3075 /* with ui and vj positive roots of the Legendre polynom */
3076 /* of degree NBPNTU and NBPNTV respectively. */
3077 /* DISOTB: Table where the terms of constraints are added */
3078 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3079 /* with ui and vj positive roots of the polynom of Legendre */
3080 /* of degree NBPNTU and NBPNTV respectively. */
3081 /* SODITB: Table where the terms of constraints are added */
3082 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3083 /* with ui and vj positive roots of the polynom of Legendre */
3084 /* of degree NBPNTU and NBPNTV respectively. */
3085 /* DIDITB: Table where the terms of constraints are added */
3086 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3087 /* with ui and vj positive roots of the polynom of Legendre */
3088 /* of degree NBPNTU and NBPNTV respectively. */
3090 /* COMMONS USED : */
3091 /* ---------------- */
3093 /* REFERENCES CALLED : */
3094 /* ----------------------- */
3096 /* DESCRIPTION/NOTES/LIMITATIONS : */
3097 /* ----------------------------------- */
3099 /* $ HISTORIQUE DES MODIFICATIONS : */
3100 /* -------------------------------- */
3101 /* 08-08-1991: RBD; Creation. */
3103 /* **********************************************************************
3106 /* Name of the routine */
3109 /* Parameter adjustments */
3111 diditb_dim1 = *nbpntu / 2 + 1;
3112 diditb_dim2 = *nbpntv / 2 + 1;
3113 diditb_offset = diditb_dim1 * diditb_dim2;
3114 diditb -= diditb_offset;
3115 disotb_dim1 = *nbpntu / 2;
3116 disotb_dim2 = *nbpntv / 2;
3117 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3118 disotb -= disotb_offset;
3119 soditb_dim1 = *nbpntu / 2;
3120 soditb_dim2 = *nbpntv / 2;
3121 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3122 soditb -= soditb_offset;
3123 sosotb_dim1 = *nbpntu / 2 + 1;
3124 sosotb_dim2 = *nbpntv / 2 + 1;
3125 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3126 sosotb -= sosotb_offset;
3127 uhermt_dim1 = (*iordru << 1) + 2;
3128 uhermt_offset = uhermt_dim1;
3129 uhermt -= uhermt_offset;
3130 fpntab_dim1 = *nbpntu;
3131 fpntab_offset = fpntab_dim1 + 1;
3132 fpntab -= fpntab_offset;
3133 ditbu2_dim1 = *nbpntv / 2 + 1;
3134 ditbu2_dim2 = *ndimen;
3135 ditbu2_offset = ditbu2_dim1 * (ditbu2_dim2 + 1);
3136 ditbu2 -= ditbu2_offset;
3137 ditbu1_dim1 = *nbpntv / 2 + 1;
3138 ditbu1_dim2 = *ndimen;
3139 ditbu1_offset = ditbu1_dim1 * (ditbu1_dim2 + 1);
3140 ditbu1 -= ditbu1_offset;
3141 sotbu2_dim1 = *nbpntv / 2 + 1;
3142 sotbu2_dim2 = *ndimen;
3143 sotbu2_offset = sotbu2_dim1 * (sotbu2_dim2 + 1);
3144 sotbu2 -= sotbu2_offset;
3145 sotbu1_dim1 = *nbpntv / 2 + 1;
3146 sotbu1_dim2 = *ndimen;
3147 sotbu1_offset = sotbu1_dim1 * (sotbu1_dim2 + 1);
3148 sotbu1 -= sotbu1_offset;
3151 ibb = AdvApp2Var_SysBase::mnfndeb_();
3153 AdvApp2Var_SysBase::mgenmsg_("MMA2CD3", 7L);
3156 /* ------------------- Discretization of polynoms of Hermit -----------
3159 ncfhu = (*iordru + 1) << 1;
3161 for (ii = 1; ii <= i__1; ++ii) {
3163 for (kk = 1; kk <= i__2; ++kk) {
3164 AdvApp2Var_MathBase::mmmpocur_(&ncfhu,
3167 &uhermt[ii * uhermt_dim1],
3169 &fpntab[kk + ii * fpntab_dim1]);
3175 /* ---- The discretizations of polynoms of constraints are subtracted ----
3178 nvroo = *nbpntv / 2;
3179 nuroo = *nbpntu / 2;
3182 for (nd = 1; nd <= i__1; ++nd) {
3184 for (ii = 1; ii <= i__2; ++ii) {
3187 for (jj = 1; jj <= i__3; ++jj) {
3188 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1];
3189 bid2 = sotbu2[jj + (nd + ii * sotbu2_dim2) * sotbu2_dim1];
3190 bid3 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1];
3191 bid4 = ditbu2[jj + (nd + ii * ditbu2_dim2) * ditbu2_dim1];
3193 for (kk = 1; kk <= i__4; ++kk) {
3194 kkp = (*nbpntu + 1) / 2 + kk;
3195 kkm = nuroo - kk + 1;
3196 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] =
3197 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
3198 - bid1 * (fpntab[kkp + ((ii << 1) - 1) *
3199 fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) *
3200 fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) *
3201 fpntab_dim1] + fpntab[kkm + (ii << 1) *
3203 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] =
3204 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
3205 - bid1 * (fpntab[kkp + ((ii << 1) - 1) *
3206 fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) *
3207 fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) *
3208 fpntab_dim1] - fpntab[kkm + (ii << 1) *
3210 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] =
3211 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
3212 - bid3 * (fpntab[kkp + ((ii << 1) - 1) *
3213 fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) *
3214 fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) *
3215 fpntab_dim1] + fpntab[kkm + (ii << 1) *
3217 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] =
3218 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
3219 - bid3 * (fpntab[kkp + ((ii << 1) - 1) *
3220 fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) *
3221 fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) *
3222 fpntab_dim1] - fpntab[kkm + (ii << 1) *
3231 /* ------------ Case when the discretization is done only on the roots */
3232 /* ---------- of Legendre polynom of uneven degree, 0 is root */
3236 if (*nbpntu % 2 == 1) {
3238 for (ii = 1; ii <= i__2; ++ii) {
3240 for (jj = 1; jj <= i__3; ++jj) {
3241 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1]
3242 * fpntab[nuroo + 1 + ((ii << 1) - 1) *
3243 fpntab_dim1] + sotbu2[jj + (nd + ii * sotbu2_dim2)
3244 * sotbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) *
3246 sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
3247 bid2 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1]
3248 * fpntab[nuroo + 1 + ((ii << 1) - 1) *
3249 fpntab_dim1] + ditbu2[jj + (nd + ii * ditbu2_dim2)
3250 * ditbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) *
3252 diditb[(jj + nd * diditb_dim2) * diditb_dim1] -= bid2;
3259 if (*nbpntv % 2 == 1) {
3261 for (ii = 1; ii <= i__2; ++ii) {
3263 for (kk = 1; kk <= i__3; ++kk) {
3264 kkp = (*nbpntu + 1) / 2 + kk;
3265 kkm = nuroo - kk + 1;
3266 bid1 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * (
3267 fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] +
3268 fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1]) +
3269 sotbu2[(nd + ii * sotbu2_dim2) * sotbu2_dim1] * (
3270 fpntab[kkp + (ii << 1) * fpntab_dim1] + fpntab[
3271 kkm + (ii << 1) * fpntab_dim1]);
3272 sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
3273 bid2 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * (
3274 fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] -
3275 fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1]) +
3276 sotbu2[(nd + ii * sotbu2_dim2) * sotbu2_dim1] * (
3277 fpntab[kkp + (ii << 1) * fpntab_dim1] - fpntab[
3278 kkm + (ii << 1) * fpntab_dim1]);
3279 diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
3286 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
3288 for (ii = 1; ii <= i__2; ++ii) {
3289 bid1 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * fpntab[
3290 nuroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbu2[(
3291 nd + ii * sotbu2_dim2) * sotbu2_dim1] * fpntab[nuroo
3292 + 1 + (ii << 1) * fpntab_dim1];
3293 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
3302 /* ------------------------------ The End -------------------------------
3307 AdvApp2Var_SysBase::mgsomsg_("MMA2CD3", 7L);
3312 //=======================================================================
3313 //function : mma2cdi_
3315 //=======================================================================
3316 int AdvApp2Var_ApproxF2var::mma2cdi_( integer *ndimen,
3342 static integer c__8 = 8;
3344 /* System generated locals */
3345 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
3346 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
3347 contr4_dim1, contr4_dim2, contr4_offset, sosotb_dim1, sosotb_dim2,
3348 sosotb_offset, diditb_dim1, diditb_dim2, diditb_offset,
3349 soditb_dim1, soditb_dim2, soditb_offset, disotb_dim1, disotb_dim2,
3352 /* Local variables */
3353 static integer ilong;
3354 static long int iofwr;
3355 static doublereal wrkar[1];
3356 static integer iszwr;
3357 static integer ibb, ier;
3358 static integer isz1, isz2, isz3, isz4;
3359 static long int ipt1, ipt2, ipt3, ipt4;
3364 /* **********************************************************************
3369 /* Discretisation on the parameters of polynomes of interpolation */
3370 /* of constraints of order IORDRE. */
3374 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3376 //* INPUT ARGUMENTS : */
3377 /* ------------------ */
3378 /* NDIMEN: Dimension of the space. */
3379 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3380 /* This is also the nb of root of Legendre polynom where discretization is done. */
3381 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
3383 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
3384 /* This is also the nb of root of Legendre polynom where discretization is done. */
3385 /* VROOTL: Table of parameters of discretisation ON (-1,1) by V.
3387 /* IORDRV: Order of constraint imposed at the extremities of iso-U */
3388 /* = 0, calculate the extremities of iso-U */
3389 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
3390 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
3391 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
3392 /* = 0, calculate the extremities of iso-V */
3393 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3394 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3395 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
3396 /* extremities of F(U0,V0) and its derivatives. */
3397 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
3398 /* extremities of F(U1,V0) and its derivatives. */
3399 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
3400 /* extremities of F(U0,V1) and its derivatives. */
3401 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
3402 /* extremities of F(U1,V1) and its derivatives. */
3403 /* SOTBU1: Table of NBPNTU/2 sums of 2 index points */
3404 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3405 /* SOTBU2: Table of NBPNTV/2 sums of 2 index points */
3406 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3407 /* DITBU1: Table of NBPNTU/2 differences of 2 index points */
3408 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3409 /* DITBU2: Table of NBPNTU/2 differences of 2 index points */
3410 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3411 /* SOTBV1: Table of NBPNTV/2 sums of 2 index points */
3412 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
3413 /* SOTBV2: Table of NBPNTV/2 sums of 2 index points */
3414 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
3415 /* DITBV1: Table of NBPNTV/2 differences of 2 index points */
3416 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
3417 /* DITBV2: Table of NBPNTV/2 differences of 2 index points */
3418 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
3419 /* SOSOTB: Preinitialized table (input/output argument). */
3420 /* DISOTB: Preinitialized table (input/output argument). */
3421 /* SODITB: Preinitialized table (input/output argument). */
3422 /* DIDITB: Preinitialized table (input/output argument) */
3424 /* ARGUMENTS DE SORTIE : */
3425 /* ------------------- */
3426 /* SOSOTB: Table where the terms of constraints are added */
3427 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3428 /* with ui and vj positive roots of the Legendre polynom */
3429 /* of degree NBPNTU and NBPNTV respectively. */
3430 /* DISOTB: Table where the terms of constraints are added */
3431 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3432 /* with ui and vj positive roots of the polynom of Legendre */
3433 /* of degree NBPNTU and NBPNTV respectively. */
3434 /* SODITB: Table where the terms of constraints are added */
3435 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3436 /* with ui and vj positive roots of the polynom of Legendre */
3437 /* of degree NBPNTU and NBPNTV respectively. */
3438 /* DIDITB: Table where the terms of constraints are added */
3439 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3440 /* with ui and vj positive roots of the polynom of Legendre */
3441 /* of degree NBPNTU and NBPNTV respectively. */
3442 /* IERCOD: = 0, OK, */
3443 /* = 1, Value or IORDRV or IORDRU is out of allowed values. */
3444 /* =13, Pb of dynamic allocation. */
3446 /* COMMONS USED : */
3447 /* ---------------- */
3449 /* REFERENCES CALLED : */
3450 /* -------------------- */
3452 /* DESCRIPTION/NOTES/LIMITATIONS : */
3453 /* ------------------------------- */
3456 /* **********************************************************************
3459 /* The name of the routine */
3462 /* Parameter adjustments */
3464 diditb_dim1 = *nbpntu / 2 + 1;
3465 diditb_dim2 = *nbpntv / 2 + 1;
3466 diditb_offset = diditb_dim1 * diditb_dim2;
3467 diditb -= diditb_offset;
3468 disotb_dim1 = *nbpntu / 2;
3469 disotb_dim2 = *nbpntv / 2;
3470 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3471 disotb -= disotb_offset;
3472 soditb_dim1 = *nbpntu / 2;
3473 soditb_dim2 = *nbpntv / 2;
3474 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3475 soditb -= soditb_offset;
3476 sosotb_dim1 = *nbpntu / 2 + 1;
3477 sosotb_dim2 = *nbpntv / 2 + 1;
3478 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3479 sosotb -= sosotb_offset;
3481 contr4_dim1 = *ndimen;
3482 contr4_dim2 = *iordru + 2;
3483 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
3484 contr4 -= contr4_offset;
3485 contr3_dim1 = *ndimen;
3486 contr3_dim2 = *iordru + 2;
3487 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
3488 contr3 -= contr3_offset;
3489 contr2_dim1 = *ndimen;
3490 contr2_dim2 = *iordru + 2;
3491 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
3492 contr2 -= contr2_offset;
3493 contr1_dim1 = *ndimen;
3494 contr1_dim2 = *iordru + 2;
3495 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
3496 contr1 -= contr1_offset;
3507 ibb = AdvApp2Var_SysBase::mnfndeb_();
3509 AdvApp2Var_SysBase::mgenmsg_("MMA2CDI", 7L);
3513 if (*iordru < -1 || *iordru > 2) {
3516 if (*iordrv < -1 || *iordrv > 2) {
3520 /* ------------------------- Set to zero --------------------------------
3523 ilong = (*nbpntu / 2 + 1) * (*nbpntv / 2 + 1) * *ndimen;
3524 AdvApp2Var_SysBase::mvriraz_(&ilong, (char *)&sosotb[sosotb_offset]);
3525 AdvApp2Var_SysBase::mvriraz_(&ilong, (char *)&diditb[diditb_offset]);
3526 ilong = *nbpntu / 2 * (*nbpntv / 2) * *ndimen;
3527 AdvApp2Var_SysBase::mvriraz_(&ilong, (char *)&soditb[soditb_offset]);
3528 AdvApp2Var_SysBase::mvriraz_(&ilong, (char *)&disotb[disotb_offset]);
3529 if (*iordru == -1 && *iordrv == -1) {
3535 isz1 = ((*iordru + 1) << 2) * (*iordru + 1);
3536 isz2 = ((*iordrv + 1) << 2) * (*iordrv + 1);
3537 isz3 = ((*iordru + 1) << 1) * *nbpntu;
3538 isz4 = ((*iordrv + 1) << 1) * *nbpntv;
3539 iszwr = isz1 + isz2 + isz3 + isz4;
3540 AdvApp2Var_SysBase::mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3549 if (*iordru >= 0 && *iordru <= 2) {
3551 /* --- Return 2*(IORDRU+1) coeff of 2*(IORDRU+1) polynoms of Hermite
3554 AdvApp2Var_ApproxF2var::mma1her_(iordru, &wrkar[ipt1], iercod);
3559 /* ---- Subract discretizations of polynoms of constraints
3562 mma2cd3_(ndimen, nbpntu, &urootl[1], nbpntv, iordru, &sotbu1[1], &
3563 sotbu2[1], &ditbu1[1], &ditbu2[1], &wrkar[ipt3], &wrkar[ipt1],
3564 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3565 disotb_offset], &diditb[diditb_offset]);
3568 if (*iordrv >= 0 && *iordrv <= 2) {
3570 /* --- Return 2*(IORDRV+1) coeff of 2*(IORDRV+1) polynoms of Hermite
3573 AdvApp2Var_ApproxF2var::mma1her_(iordrv, &wrkar[ipt2], iercod);
3578 /* ---- Subtract discretisations of polynoms of constraint
3581 mma2cd2_(ndimen, nbpntu, nbpntv, &vrootl[1], iordrv, &sotbv1[1], &
3582 sotbv2[1], &ditbv1[1], &ditbv2[1], &wrkar[ipt4], &wrkar[ipt2],
3583 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3584 disotb_offset], &diditb[diditb_offset]);
3587 /* --------------- Subtract constraints of corners ----------------
3590 if (*iordru >= 0 && *iordrv >= 0) {
3591 mma2cd1_(ndimen, nbpntu, &urootl[1], nbpntv, &vrootl[1], iordru,
3592 iordrv, &contr1[contr1_offset], &contr2[contr2_offset], &
3593 contr3[contr3_offset], &contr4[contr4_offset], &wrkar[ipt3], &
3594 wrkar[ipt4], &wrkar[ipt1], &wrkar[ipt2], &sosotb[
3595 sosotb_offset], &soditb[soditb_offset], &disotb[disotb_offset]
3596 , &diditb[diditb_offset]);
3600 /* ------------------------------ The End -------------------------------
3602 /* --> IORDRE is not within the autorised diapason. */
3606 /* --> PB of dynamic allocation. */
3613 AdvApp2Var_SysBase::mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3618 AdvApp2Var_SysBase::maermsg_("MMA2CDI", iercod, 7L);
3620 AdvApp2Var_SysBase::mgsomsg_("MMA2CDI", 7L);
3625 //=======================================================================
3626 //function : mma2ce1_
3628 //=======================================================================
3629 int AdvApp2Var_ApproxF2var::mma2ce1_(integer *numdec,
3657 static integer c__8 = 8;
3659 /* System generated locals */
3660 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3661 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3662 diditb_dim1, diditb_dim2, diditb_offset, patjac_dim1, patjac_dim2,
3665 /* Local variables */
3666 static logical ldbg;
3667 static long int iofwr;
3668 static doublereal wrkar[1];
3669 static integer iszwr;
3671 static integer isz1, isz2, isz3, isz4, isz5, isz6, isz7;
3672 static long int ipt1, ipt2, ipt3, ipt4, ipt5, ipt6, ipt7;
3676 /* **********************************************************************
3681 /* Calculation of coefficients of polynomial approximation of degree */
3682 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3683 /* discretization on roots of Legendre polynom of degree */
3684 /* NBPNTU by U and NBPNTV by V. */
3688 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&POLYNOME,&ERREUR */
3690 /* INPUT ARGUMENTS : */
3691 /* ------------------ */
3692 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3693 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3694 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3695 /* directions simultaneously (cutting by V is preferable). */
3696 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3697 /* directions simultaneously (cutting by U is preferable). */
3698 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3699 /* of cutting Vj). */
3700 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3701 /* of cutting Ui). */
3702 /* = 0, It is not POSSIBLE to cut anything */
3703 /* NDIMEN: Dimension of the space. */
3704 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3705 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3706 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3707 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3708 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3709 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3710 /* NDJACU: Max degree of the polynom of approximation by U. */
3711 /* The representation in the orthogonal base starts from degree */
3712 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3713 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3714 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3715 /* NDJACV: Max degree of the polynom of approximation by V. */
3716 /* The representation in the orthogonal base starts from degree */
3717 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3718 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3719 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3720 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3721 /* to the step of constraints C0, C1 or C2. */
3722 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3723 /* to the step of constraints C0, C1 or C2. */
3724 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3725 /* calculated the coefficients of integration by u */
3726 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3727 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3728 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3729 /* calculated the coefficients of integration by u */
3730 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3731 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3732 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3733 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3734 /* with ui and vj - positive roots of the Legendre polynom */
3735 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3736 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3737 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3738 /* SOSOTB(0,0) contains F(0,0). */
3739 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3740 /* with ui and vj positive roots of Legendre polynom */
3741 /* of degree NBPNTU and NBPNTV respectively. */
3742 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3743 /* with ui and vj positive roots of Legendre polynom */
3744 /* of degree NBPNTU and NBPNTV respectively. */
3745 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3746 /* with ui and vj positive roots of Legendre polynom */
3747 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3748 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3749 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3751 /* OUTPUT ARGUMENTS */
3752 /* --------------- */
3753 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
3754 /* of F(u,v) with eventually taking into account of */
3755 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
3756 /* This table contains other coeff if ITYDEC = 0. */
3757 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
3758 /* on each of sub-spaces SI ITYDEC = 0. */
3759 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
3760 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
3761 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
3762 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
3763 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
3764 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
3765 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
3766 /* = 3, it is NECESSARY to cut both by U AND by V. */
3767 /* IERCOD: Error code. */
3768 /* = 0, Everything is OK. */
3769 /* = -1, There is the best possible solution, but the */
3770 /* user tolerance is not satisfactory (3*only) */
3771 /* = 1, Incoherent entries. */
3773 /* COMMONS USED : */
3774 /* ---------------- */
3776 /* REFERENCES CALLED : */
3777 /* --------------------- */
3779 /* DESCRIPTION/NOTES/LIMITATIONS : */
3780 /* ------------------------------- */
3783 /* **********************************************************************
3785 /* Name of the routine */
3788 /* --------------------------- Initialisations --------------------------
3791 /* Parameter adjustments */
3796 patjac_dim1 = *ndjacu + 1;
3797 patjac_dim2 = *ndjacv + 1;
3798 patjac_offset = patjac_dim1 * patjac_dim2;
3799 patjac -= patjac_offset;
3800 diditb_dim1 = *nbpntu / 2 + 1;
3801 diditb_dim2 = *nbpntv / 2 + 1;
3802 diditb_offset = diditb_dim1 * diditb_dim2;
3803 diditb -= diditb_offset;
3804 soditb_dim1 = *nbpntu / 2;
3805 soditb_dim2 = *nbpntv / 2;
3806 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3807 soditb -= soditb_offset;
3808 disotb_dim1 = *nbpntu / 2;
3809 disotb_dim2 = *nbpntv / 2;
3810 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3811 disotb -= disotb_offset;
3812 sosotb_dim1 = *nbpntu / 2 + 1;
3813 sosotb_dim2 = *nbpntv / 2 + 1;
3814 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3815 sosotb -= sosotb_offset;
3818 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
3820 AdvApp2Var_SysBase::mgenmsg_("MMA2CE1", 7L);
3825 isz1 = (*nbpntu / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1);
3826 isz2 = (*nbpntv / 2 + 1) * (*ndjacv - ((*iordrv + 1) << 1) + 1);
3827 isz3 = (*nbpntv / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3828 isz4 = *nbpntv / 2 * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3829 isz5 = *ndjacu + 1 - ((*iordru + 1) << 1);
3830 isz6 = *ndjacv + 1 - ((*iordrv + 1) << 1);
3831 isz7 = *ndimen << 2;
3832 iszwr = isz1 + isz2 + isz3 + isz4 + isz5 + isz6 + isz7;
3833 AdvApp2Var_SysBase::mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3845 /* ----------------- Return Gauss coefficients of integration ----------------
3848 AdvApp2Var_ApproxF2var::mmapptt_(ndjacu, nbpntu, iordru, &wrkar[ipt1], iercod);
3852 AdvApp2Var_ApproxF2var::mmapptt_(ndjacv, nbpntv, iordrv, &wrkar[ipt2], iercod);
3857 /* ------------------- Return max polynoms of Jacobi ------------
3860 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacu, iordru, &wrkar[ipt5]);
3861 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacv, iordrv, &wrkar[ipt6]);
3863 /* ------ Calculate the coefficients and their contribution to the error ----
3866 mma2ce2_(numdec, ndimen, nbsesp, &ndimse[1], ndminu, ndminv, ndguli,
3867 ndgvli, ndjacu, ndjacv, iordru, iordrv, nbpntu, nbpntv, &epsapr[1]
3868 , &sosotb[sosotb_offset], &disotb[disotb_offset], &soditb[
3869 soditb_offset], &diditb[diditb_offset], &wrkar[ipt1], &wrkar[ipt2]
3870 , &wrkar[ipt5], &wrkar[ipt6], &wrkar[ipt7], &wrkar[ipt3], &wrkar[
3871 ipt4], &patjac[patjac_offset], &errmax[1], &errmoy[1], ndegpu,
3872 ndegpv, itydec, iercod);
3878 /* ------------------------------ The end -------------------------------
3887 AdvApp2Var_SysBase::mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3892 AdvApp2Var_SysBase::maermsg_("MMA2CE1", iercod, 7L);
3894 AdvApp2Var_SysBase::mgsomsg_("MMA2CE1", 7L);
3899 //=======================================================================
3900 //function : mma2ce2_
3902 //=======================================================================
3903 int mma2ce2_(integer *numdec,
3938 /* System generated locals */
3939 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3940 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3941 diditb_dim1, diditb_dim2, diditb_offset, gssutb_dim1, gssvtb_dim1,
3942 chpair_dim1, chpair_dim2, chpair_offset, chimpr_dim1,
3943 chimpr_dim2, chimpr_offset, patjac_dim1, patjac_dim2,
3944 patjac_offset, vecerr_dim1, vecerr_offset, i__1, i__2, i__3, i__4;
3946 /* Local variables */
3947 static logical ldbg;
3948 static integer idim, igsu, minu, minv, maxu, maxv, igsv;
3949 static doublereal vaux[3];
3950 static integer i2rdu, i2rdv, ndses, nd, ii, jj, kk, nu, nv;
3951 static doublereal zu, zv;
3952 static integer nu1, nv1;
3954 /* **********************************************************************
3958 /* Calculation of coefficients of polynomial approximation of degree */
3959 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3960 /* discretization on roots of Legendre polynom of degree */
3961 /* NBPNTU by U and NBPNTV by V. */
3965 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&COEFFICIENT,&POLYNOME */
3967 /* INPUT ARGUMENTS : */
3968 /* ------------------ */
3969 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3970 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3971 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3972 /* directions simultaneously (cutting by V is preferable). */
3973 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3974 /* directions simultaneously (cutting by U is preferable). */
3975 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3976 /* of cutting Vj). */
3977 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3978 /* of cutting Ui). */
3979 /* = 0, It is not POSSIBLE to cut anything */
3980 /* NDIMEN: Total dimension of the space. */
3981 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3982 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3983 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3984 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3985 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3986 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3987 /* NDJACU: Max degree of the polynom of approximation by U. */
3988 /* The representation in the orthogonal base starts from degree */
3989 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3990 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3991 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3992 /* NDJACV: Max degree of the polynom of approximation by V. */
3993 /* The representation in the orthogonal base starts from degree */
3994 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3995 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3996 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3997 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3998 /* to the step of constraints C0, C1 or C2. */
3999 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
4000 /* to the step of constraints C0, C1 or C2. */
4001 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
4002 /* calculated the coefficients of integration by u */
4003 /* by Gauss method. It is required that NBPNTU = 30, 40, */
4004 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
4005 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
4006 /* calculated the coefficients of integration by u */
4007 /* by Gauss method. It is required that NBPNTV = 30, 40, */
4008 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
4009 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
4010 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
4011 /* with ui and vj - positive roots of the Legendre polynom */
4012 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
4013 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
4014 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
4015 /* SOSOTB(0,0) contains F(0,0). */
4016 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
4017 /* with ui and vj positive roots of Legendre polynom */
4018 /* of degree NBPNTU and NBPNTV respectively. */
4019 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
4020 /* with ui and vj positive roots of Legendre polynom */
4021 /* of degree NBPNTU and NBPNTV respectively. */
4022 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
4023 /* with ui and vj positive roots of Legendre polynom */
4024 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
4025 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
4026 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
4027 /* GSSUTB: Table of coefficients of integration by Gauss method */
4028 /* by U: i varies from 0 to NBPNTU/2 and k varies from 0 to */
4029 /* NDJACU-2*(IORDRU+1). */
4030 /* GSSVTB: Table of coefficients of integration by Gauss method */
4031 /* by V: i varies from 0 to NBPNTV/2 and k varies from 0 to */
4032 /* NDJACV-2*(IORDRV+1). */
4033 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
4034 /* from degree 0 to degree NDJACU - 2*(IORDRU+1) */
4035 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
4036 /* from degree 0 to degree NDJACV - 2*(IORDRV+1) */
4038 /* OUTPUT ARGUMENTS : */
4039 /* ------------------- */
4040 /* VECERR: Auxiliary table. */
4041 /* CHPAIR: Auxiliary table of terms connected to degree NDJACU by U */
4042 /* to calculate the coeff. of approximation of EVEN degree by V. */
4043 /* CHIMPR: Auxiliary table of terms connected to degree NDJACU by U */
4044 /* to calculate the coeff. of approximation of UNEVEN degree by V. */
4045 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
4046 /* of F(u,v) with eventually taking into account of */
4047 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
4048 /* This table contains other coeff if ITYDEC = 0. */
4049 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
4050 /* on each of sub-spaces SI ITYDEC = 0. */
4051 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
4052 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
4053 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
4054 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
4055 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
4056 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
4057 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
4058 /* = 3, it is NECESSARY to cut both by U AND by V. */
4059 /* IERCOD: Error code. */
4060 /* = 0, Everything is OK. */
4061 /* = -1, There is the best possible solution, but the */
4062 /* user tolerance is not satisfactory (3*only) */
4063 /* = 1, Incoherent entries. */
4065 /* COMMONS USED : */
4066 /* ---------------- */
4068 /* REFERENCES CALLED : */
4069 /* --------------------- */
4071 /* DESCRIPTION/NOTES/LIMITATIONS : */
4073 /* **********************************************************************
4075 /* Name of the routine */
4078 /* --------------------------- Initialisations --------------------------
4081 /* Parameter adjustments */
4082 vecerr_dim1 = *ndimen;
4083 vecerr_offset = vecerr_dim1 + 1;
4084 vecerr -= vecerr_offset;
4089 patjac_dim1 = *ndjacu + 1;
4090 patjac_dim2 = *ndjacv + 1;
4091 patjac_offset = patjac_dim1 * patjac_dim2;
4092 patjac -= patjac_offset;
4093 gssutb_dim1 = *nbpntu / 2 + 1;
4094 chimpr_dim1 = *nbpntv / 2;
4095 chimpr_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4096 chimpr_offset = chimpr_dim1 * chimpr_dim2 + 1;
4097 chimpr -= chimpr_offset;
4098 chpair_dim1 = *nbpntv / 2 + 1;
4099 chpair_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4100 chpair_offset = chpair_dim1 * chpair_dim2;
4101 chpair -= chpair_offset;
4102 gssvtb_dim1 = *nbpntv / 2 + 1;
4103 diditb_dim1 = *nbpntu / 2 + 1;
4104 diditb_dim2 = *nbpntv / 2 + 1;
4105 diditb_offset = diditb_dim1 * diditb_dim2;
4106 diditb -= diditb_offset;
4107 soditb_dim1 = *nbpntu / 2;
4108 soditb_dim2 = *nbpntv / 2;
4109 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
4110 soditb -= soditb_offset;
4111 disotb_dim1 = *nbpntu / 2;
4112 disotb_dim2 = *nbpntv / 2;
4113 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
4114 disotb -= disotb_offset;
4115 sosotb_dim1 = *nbpntu / 2 + 1;
4116 sosotb_dim2 = *nbpntv / 2 + 1;
4117 sosotb_offset = sosotb_dim1 * sosotb_dim2;
4118 sosotb -= sosotb_offset;
4121 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4123 AdvApp2Var_SysBase::mgenmsg_("MMA2CE2", 7L);
4125 /* --> A priori everything is OK */
4127 /* --> test of inputs */
4128 if (*numdec < 0 || *numdec > 5) {
4131 if ((*iordru << 1) + 1 > *ndminu) {
4134 if (*ndminu > *ndguli) {
4137 if (*ndguli >= *ndjacu) {
4140 if ((*iordrv << 1) + 1 > *ndminv) {
4143 if (*ndminv > *ndgvli) {
4146 if (*ndgvli >= *ndjacv) {
4149 /* --> A priori, no cuts to be done */
4151 /* --> Min. degrees to return: NDMINU,NDMINV */
4154 /* --> For the moment, max errors are null */
4155 AdvApp2Var_SysBase::mvriraz_(nbsesp, (char *)&errmax[1]);
4157 AdvApp2Var_SysBase::mvriraz_(&nd, (char *)&vecerr[vecerr_offset]);
4158 /* --> and the square, too. */
4159 nd = (*ndjacu + 1) * (*ndjacv + 1) * *ndimen;
4160 AdvApp2Var_SysBase::mvriraz_(&nd, (char *)&patjac[patjac_offset]);
4162 i2rdu = (*iordru + 1) << 1;
4163 i2rdv = (*iordrv + 1) << 1;
4165 /* **********************************************************************
4167 /* -------------------- HERE IT IS POSSIBLE TO CUT ----------------------
4169 /* **********************************************************************
4172 if (*numdec > 0 && *numdec <= 5) {
4174 /* ******************************************************************
4176 /* ---------------------- Calculate coeff of zone 4 -------------
4190 /* ---------------- Calculate the terms connected to degree by U ---------
4194 for (nd = 1; nd <= i__1; ++nd) {
4196 for (kk = minu; kk <= i__2; ++kk) {
4198 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4199 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4200 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4201 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4202 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4203 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4204 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4210 /* ------------------- Calculate the coefficients of PATJAC ------------
4213 igsu = minu - i2rdu;
4215 for (jj = minv; jj <= i__1; ++jj) {
4218 for (nd = 1; nd <= i__2; ++nd) {
4219 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4220 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4221 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4222 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4223 patjac_dim2) * patjac_dim1]);
4227 /* ----- Contribution of calculated terms to the approximation error */
4228 /* for terms (I,J) with MINU <= I <= MAXU, J fixe. */
4232 for (nd = 1; nd <= i__2; ++nd) {
4234 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &jj, &jj,
4235 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4236 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4237 &vecerr[nd + (vecerr_dim1 << 2)]);
4238 if (vecerr[nd + (vecerr_dim1 << 2)] > epsapr[nd]) {
4247 /* ******************************************************************
4249 /* ---------------------- Calculate the coeff of zone 2 -------------
4252 minu = (*iordru + 1) << 1;
4257 /* --> If zone 2 is empty, pass to zone 3. */
4258 /* VECERR(ND,2) was already set to zero. */
4263 /* ---------------- Calculate the terms connected to degree by U ------------
4267 for (nd = 1; nd <= i__1; ++nd) {
4269 for (kk = minu; kk <= i__2; ++kk) {
4271 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4272 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4273 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4274 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4275 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4276 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4277 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4283 /* ------------------- Calculate the coefficients of PATJAC ------------
4286 igsu = minu - i2rdu;
4288 for (jj = minv; jj <= i__1; ++jj) {
4291 for (nd = 1; nd <= i__2; ++nd) {
4292 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4293 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4294 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4295 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4296 patjac_dim2) * patjac_dim1]);
4302 /* -----Contribution of calculated terms to the approximation error */
4303 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4307 for (nd = 1; nd <= i__1; ++nd) {
4309 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4310 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4311 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4312 vecerr[nd + (vecerr_dim1 << 1)]);
4317 /* ******************************************************************
4319 /* ---------------------- Calculation of coeff of zone 3 -------------
4325 minv = (*iordrv + 1) << 1;
4328 /* -> If zone 3 is empty, pass to the test of cutting. */
4329 /* VECERR(ND,3) was already set to zero */
4334 /* ----------- The terms connected to the degree by U are already calculated -----
4336 /* ------------------- Calculation of coefficients of PATJAC ------------
4339 igsu = minu - i2rdu;
4341 for (jj = minv; jj <= i__1; ++jj) {
4344 for (nd = 1; nd <= i__2; ++nd) {
4345 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4346 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4347 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4348 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4349 patjac_dim2) * patjac_dim1]);
4355 /* ----- Contribution of calculated terms to the approximation error
4356 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV. */
4360 for (nd = 1; nd <= i__1; ++nd) {
4362 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4363 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4364 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4365 vecerr[nd + vecerr_dim1 * 3]);
4370 /* ******************************************************************
4372 /* --------------------------- Tests of cutting ---------------------
4377 for (nd = 1; nd <= i__1; ++nd) {
4378 vaux[0] = vecerr[nd + (vecerr_dim1 << 1)];
4379 vaux[1] = vecerr[nd + (vecerr_dim1 << 2)];
4380 vaux[2] = vecerr[nd + vecerr_dim1 * 3];
4382 errmax[nd] = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4383 if (errmax[nd] > epsapr[nd]) {
4385 zv = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4386 zu = AdvApp2Var_MathBase::mzsnorm_(&ii, &vaux[1]);
4387 if (zu > epsapr[nd] && zv > epsapr[nd]) {
4399 /* ******************************************************************
4401 /* --- OK, the square is valid, the coeff of zone 1 are calculated
4404 minu = (*iordru + 1) << 1;
4406 minv = (*iordrv + 1) << 1;
4409 /* --> If zone 1 is empty, pass to the calculation of Max and Average error. */
4410 if (minu > maxu || minv > maxv) {
4414 /* ----------- The terms connected to degree by U are already calculated -----
4416 /* ------------------- Calculate the coefficients of PATJAC ------------
4419 igsu = minu - i2rdu;
4421 for (jj = minv; jj <= i__1; ++jj) {
4424 for (nd = 1; nd <= i__2; ++nd) {
4425 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4426 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4427 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4428 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4429 patjac_dim2) * patjac_dim1]);
4435 /* --------------- Now the degree is maximally lowered --------
4440 i__1 = 1, i__2 = (*iordru << 1) + 1, i__1 = max(i__1,i__2);
4441 minu = max(i__1,*ndminu);
4444 i__1 = 1, i__2 = (*iordrv << 1) + 1, i__1 = max(i__1,i__2);
4445 minv = max(i__1,*ndminv);
4449 for (nd = 1; nd <= i__1; ++nd) {
4451 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4452 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4453 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4454 patjac_dim2 * patjac_dim1], &epsapr[nd], &vecerr[
4455 vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4463 /* --> Calculate the average error. */
4464 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4465 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4468 /* --> Set to 0.D0 the rejected coeffs. */
4469 i__2 = idim + ndses - 1;
4470 for (ii = idim; ii <= i__2; ++ii) {
4472 for (jj = nv1; jj <= i__3; ++jj) {
4474 for (kk = nu1; kk <= i__4; ++kk) {
4475 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4484 /* --> Return the nb of coeffs of approximation. */
4485 *ndegpu = max(*ndegpu,nu);
4486 *ndegpv = max(*ndegpv,nv);
4491 /* ******************************************************************
4493 /* -------------------- IT IS NOT POSSIBLE TO CUT -------------------
4495 /* ******************************************************************
4499 minu = (*iordru + 1) << 1;
4501 minv = (*iordrv + 1) << 1;
4504 /* ---------------- Calculate the terms connected to the degree by U ------------
4508 for (nd = 1; nd <= i__1; ++nd) {
4510 for (kk = minu; kk <= i__2; ++kk) {
4512 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4513 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4514 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4515 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4516 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4517 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4518 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4522 /* ---------------------- Calculate all coefficients -------
4525 igsu = minu - i2rdu;
4527 for (jj = minv; jj <= i__2; ++jj) {
4529 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4530 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4531 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4532 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4533 patjac_dim2) * patjac_dim1]);
4539 /* ----- Contribution of calculated terms to the approximation error
4540 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4544 for (nd = 1; nd <= i__1; ++nd) {
4546 minu = (*iordru + 1) << 1;
4550 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4551 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4552 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4556 minv = (*iordrv + 1) << 1;
4559 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4560 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4561 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4565 /* ---------------------------- IF ERRMAX > EPSAPR, stop --------
4568 if (errmax[nd] > epsapr[nd]) {
4573 /* ------------- Otherwise, try to remove again the coeff
4578 i__2 = 1, i__3 = (*iordru << 1) + 1, i__2 = max(i__2,i__3);
4579 minu = max(i__2,*ndminu);
4582 i__2 = 1, i__3 = (*iordrv << 1) + 1, i__2 = max(i__2,i__3);
4583 minv = max(i__2,*ndminv);
4585 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4586 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &
4587 maxv, iordru, iordrv, xmaxju, xmaxjv, &patjac[
4588 idim * patjac_dim2 * patjac_dim1], &epsapr[nd], &
4589 vecerr[vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4596 /* --------------------- Calculate the average error -------------
4601 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4602 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4605 /* --------------------- Set to 0.D0 the rejected coeffs ----------
4608 i__2 = idim + ndses - 1;
4609 for (ii = idim; ii <= i__2; ++ii) {
4611 for (jj = nv1; jj <= i__3; ++jj) {
4613 for (kk = nu1; kk <= i__4; ++kk) {
4614 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4623 /* --------------- Return the nb of coeff of approximation ---
4626 *ndegpu = max(*ndegpu,nu);
4627 *ndegpv = max(*ndegpv,nv);
4635 /* ------------------------------ The end -------------------------------
4637 /* --> Error in inputs */
4642 /* --------- Management of cuts, it is required 0 < NUMDEC <= 5 -------
4645 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4647 if (*numdec <= 0 || *numdec > 5) {
4656 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4658 if (*numdec <= 0 || *numdec > 5) {
4667 /* --> Here it is possible and necessary to cut, choose by 4 if it is possible */
4669 if (*numdec <= 0 || *numdec > 5) {
4674 } else if (*numdec == 2 || *numdec == 4) {
4676 } else if (*numdec == 1 || *numdec == 3) {
4684 AdvApp2Var_SysBase::maermsg_("MMA2CE2", iercod, 7L);
4686 AdvApp2Var_SysBase::mgsomsg_("MMA2CE2", 7L);
4691 //=======================================================================
4692 //function : mma2cfu_
4694 //=======================================================================
4695 int mma2cfu_(integer *ndujac,
4707 /* System generated locals */
4708 integer sosotb_dim1, disotb_dim1, disotb_offset, soditb_dim1,
4709 soditb_offset, diditb_dim1, i__1, i__2;
4711 /* Local variables */
4712 static logical ldbg;
4713 static integer nptu2, nptv2, ii, jj;
4714 static doublereal bid0, bid1, bid2;
4717 /* **********************************************************************
4722 /* Calculate the terms connected to degree NDUJAC by U of the polynomial approximation */
4723 /* of function F(u,v), starting from its discretisation
4724 /* on the roots of Legendre polynom of degree */
4725 /* NBPNTU by U and NBPNTV by V. */
4729 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4731 /* INPUT ARGUMENTSE : */
4732 /* ------------------ */
4733 /* NDUJAC: Fixed degree by U for which the terms */
4734 /* allowing to obtain the Legendre or Jacobi coeff*/
4735 /* of even or uneven degree by V are calculated. */
4736 /* NBPNTU: Degree of Legendre polynom on the roots which of */
4737 /* the coefficients of integration by U are calculated */
4738 /* by Gauss method. It is required that NBPNTU = 30, 40, 50 or 61. */
4739 /* NBPNTV: Degree of Legendre polynom on the roots which of */
4740 /* the coefficients of integration by V are calculated */
4741 /* by Gauss method. It is required that NBPNTV = 30, 40, 50 or 61. */
4742 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
4743 /* with ui and vj positive roots of Legendre polynom */
4744 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4745 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
4746 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
4747 /* SOSOTB(0,0) contains F(0,0). */
4748 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
4749 /* with ui and vj positive roots of Legendre polynom */
4750 /* of degree NBPNTU and NBPNTV respectively. */
4751 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
4752 /* with ui and vj positive roots of Legendre polynom */
4753 /* of degree NBPNTU and NBPNTV respectively. */
4754 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
4755 /* avec ui and vj positive roots of Legendre polynom */
4756 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4757 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
4758 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
4759 /* GSSUTB: Table of coefficients of integration by Gauss method */
4760 /* Gauss by U for fixed NDUJAC : i varies from 0 to NBPNTU/2. */
4762 /* OUTPUT ARGUMENTS : */
4763 /* ------------------- */
4764 /* CHPAIR: Table of terms connected to degree NDUJAC by U to calculate the */
4765 /* coeff. of the approximation of EVEN degree by V. */
4766 /* CHIMPR: Table of terms connected to degree NDUJAC by U to calculate */
4767 /* the coeff. of approximation of UNEVEN degree by V. */
4769 /* COMMONS USED : */
4770 /* ---------------- */
4772 /* REFERENCES CALLED : */
4773 /* ----------------------- */
4775 /* DESCRIPTION/NOTES/LIMITATIONS : */
4776 /* ----------------------------------- */
4780 /* **********************************************************************
4782 /* Name of the routine */
4785 /* --------------------------- Initialisations --------------------------
4788 /* Parameter adjustments */
4790 diditb_dim1 = *nbpntu / 2 + 1;
4791 soditb_dim1 = *nbpntu / 2;
4792 soditb_offset = soditb_dim1 + 1;
4793 soditb -= soditb_offset;
4794 disotb_dim1 = *nbpntu / 2;
4795 disotb_offset = disotb_dim1 + 1;
4796 disotb -= disotb_offset;
4797 sosotb_dim1 = *nbpntu / 2 + 1;
4800 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4802 AdvApp2Var_SysBase::mgenmsg_("MMA2CFU", 7L);
4805 nptu2 = *nbpntu / 2;
4806 nptv2 = *nbpntv / 2;
4808 /* **********************************************************************
4810 /* CALCULATE COEFFICIENTS BY U */
4812 /* ----------------- Calculate coefficients of even degree --------------
4815 if (*ndujac % 2 == 0) {
4817 for (jj = 1; jj <= i__1; ++jj) {
4821 for (ii = 1; ii <= i__2; ++ii) {
4823 bid1 += sosotb[ii + jj * sosotb_dim1] * bid0;
4824 bid2 += soditb[ii + jj * soditb_dim1] * bid0;
4832 /* --------------- Calculate coefficients of uneven degree ----------
4837 for (jj = 1; jj <= i__1; ++jj) {
4841 for (ii = 1; ii <= i__2; ++ii) {
4843 bid1 += disotb[ii + jj * disotb_dim1] * bid0;
4844 bid2 += diditb[ii + jj * diditb_dim1] * bid0;
4853 /* ------- Add terms connected to the supplementary root (0.D0) ------
4854 /* ----------- of Legendre polynom of uneven degree NBPNTU -----------
4856 /* --> Only even NDUJAC terms are modified as GSSUTB(0) = 0 */
4857 /* when NDUJAC is uneven. */
4859 if (*nbpntu % 2 != 0 && *ndujac % 2 == 0) {
4862 for (jj = 1; jj <= i__1; ++jj) {
4863 chpair[jj] += sosotb[jj * sosotb_dim1] * bid0;
4864 chimpr[jj] += diditb[jj * diditb_dim1] * bid0;
4869 /* ------ Calculate the terms connected to supplementary roots (0.D0) ------
4871 /* ----------- of Legendre polynom of uneven degree NBPNTV -----------
4874 if (*nbpntv % 2 != 0) {
4875 /* --> Only CHPAIR terms are calculated as GSSVTB(0,IH-IDEBV)=0
4877 /* when IH is uneven (see MMA2CFV). */
4879 if (*ndujac % 2 == 0) {
4882 for (ii = 1; ii <= i__1; ++ii) {
4883 bid1 += sosotb[ii] * gssutb[ii];
4890 for (ii = 1; ii <= i__1; ++ii) {
4891 bid1 += diditb[ii] * gssutb[ii];
4896 if (*nbpntu % 2 != 0) {
4897 chpair[0] += sosotb[0] * gssutb[0];
4901 /* ------------------------------ The end -------------------------------
4905 AdvApp2Var_SysBase::mgsomsg_("MMA2CFU", 7L);
4910 //=======================================================================
4911 //function : mma2cfv_
4913 //=======================================================================
4914 int mma2cfv_(integer *ndvjac,
4924 /* System generated locals */
4925 integer chpair_dim1, chpair_offset, chimpr_dim1, chimpr_offset,
4926 patjac_offset, i__1, i__2;
4928 /* Local variables */
4929 static logical ldbg;
4930 static integer nptv2, ii, jj;
4931 static doublereal bid1;
4934 /* **********************************************************************
4939 /* Calculate the coefficients of polynomial approximation of F(u,v)
4940 /* of degree NDVJAC by V and of degree by U varying from MINDGU to MAXDGU.
4945 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4947 /* INPUT ARGUMENTS : */
4948 /* ------------------ */
4950 /* NDVJAC: Degree of the polynom of approximation by V. */
4951 /* The representation in the orthogonal base starts from degre 0.
4952 /* The polynomial base is the base of Jacobi of order -1 */
4953 /* (Legendre), 0, 1 or 2 */
4954 /* MINDGU: Degree minimum by U of coeff. to calculate. */
4955 /* MAXDGU: Degree maximum by U of coeff. to calculate. */
4956 /* NBPNTV: Degree of the Legendre polynom on the roots which of */
4957 /* the coefficients of integration by V are calculated */
4958 /* by Gauss method. It is reqired that NBPNTV = 30, 40, 50 or 61 and NDVJAC < NBPNTV. */
4959 /* GSSVTB: Table of coefficients of integration by Gauss method */
4960 /* by V for NDVJAC fixed: j varies from 0 to NBPNTV/2. */
4961 /* CHPAIR: Table of terms connected to degrees from MINDGU to MAXDGU by U to
4962 /* calculate the coeff. of approximation of EVEN degree NDVJAC by V. */
4963 /* CHIMPR: Table of terms connected to degrees from MINDGU to MAXDGU by U to
4964 /* calculate the coeff. of approximation of UNEVEN degree NDVJAC by V. */
4966 /* OUTPUT ARGUMENTS : */
4967 /* ------------------- */
4968 /* PATJAC: Table of coefficients by U of the polynom of approximation */
4969 /* P(u,v) of degree MINDGU to MAXDGU by U and NDVJAC by V. */
4971 /* COMMONS USED : */
4972 /* -------------- */
4974 /* REFERENCES CALLED : */
4975 /* --------------------- */
4977 /* DESCRIPTION/NOTES/LIMITATIONS : */
4978 /* ------------------------------- */
4980 /* **********************************************************************
4982 /* Name of the routine */
4985 /* --------------------------- Initialisations --------------------------
4988 /* Parameter adjustments */
4989 patjac_offset = *mindgu;
4990 patjac -= patjac_offset;
4991 chimpr_dim1 = *nbpntv / 2;
4992 chimpr_offset = chimpr_dim1 * *mindgu + 1;
4993 chimpr -= chimpr_offset;
4994 chpair_dim1 = *nbpntv / 2 + 1;
4995 chpair_offset = chpair_dim1 * *mindgu;
4996 chpair -= chpair_offset;
4999 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5001 AdvApp2Var_SysBase::mgenmsg_("MMA2CFV", 7L);
5003 nptv2 = *nbpntv / 2;
5005 /* --------- Calculate the coefficients for even degree NDVJAC ----------
5008 if (*ndvjac % 2 == 0) {
5010 for (ii = *mindgu; ii <= i__1; ++ii) {
5013 for (jj = 1; jj <= i__2; ++jj) {
5014 bid1 += chpair[jj + ii * chpair_dim1] * gssvtb[jj];
5021 /* -------- Calculate the coefficients for uneven degree NDVJAC -----
5026 for (ii = *mindgu; ii <= i__1; ++ii) {
5029 for (jj = 1; jj <= i__2; ++jj) {
5030 bid1 += chimpr[jj + ii * chimpr_dim1] * gssvtb[jj];
5038 /* ------- Add terms connected to the supplementary root (0.D0) ----- */
5039 /* --------of the Legendre polynom of uneven degree NBPNTV --------- */
5041 if (*nbpntv % 2 != 0 && *ndvjac % 2 == 0) {
5044 for (ii = *mindgu; ii <= i__1; ++ii) {
5045 patjac[ii] += bid1 * chpair[ii * chpair_dim1];
5050 /* ------------------------------ The end -------------------------------
5054 AdvApp2Var_SysBase::mgsomsg_("MMA2CFV", 7L);
5059 //=======================================================================
5060 //function : mma2ds1_
5062 //=======================================================================
5063 int AdvApp2Var_ApproxF2var::mma2ds1_(integer *ndimen,
5093 /* System generated locals */
5094 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5095 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5096 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5097 fpntab_offset, i__1;
5099 /* Local variables */
5100 static logical ldbg;
5101 static integer ibid1, ibid2, iuouv, nd;
5102 static integer isz1, isz2;
5106 /* **********************************************************************
5111 /* Discretisation of function F(u,v) on the roots of Legendre polynoms. */
5115 /* FONCTION&,DISCRETISATION,&POINT */
5117 /* INPUT ARGUMENTS : */
5118 /* ------------------ */
5119 /* NDIMEN: Dimension of the space. */
5120 /* UINTFN: Limits of the interval of definition by u of the function */
5121 /* to be processed: (UINTFN(1),UINTFN(2)). */
5122 /* VINTFN: Limits of the interval of definition by v of the function */
5123 /* to be processed: (VINTFN(1),VINTFN(2)). */
5124 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5125 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5126 /* FONCNP is discretized by u. */
5127 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5128 /* FONCNP is discretized by v. */
5129 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5130 /* of Legendre of degree NBPNTU defined on (-1,1). */
5131 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5132 /* of Legendre of degree NBPNTV defined on (-1,1). */
5133 /* ISOFAV: Shows the type of iso of F(u,v) to be extracted to improve */
5134 /* the rapidity of calculation (has no influence on the form */
5136 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5137 /* with fixed u (with NBPNTV values different from v). */
5138 /* = 2, shows that it is necessaty to calculate the points of F(u,v) */
5139 /* with fixed v (with NBPNTU values different from u). */
5140 /* SOSOTB: Preinitialized table (input/output argument). */
5141 /* DISOTB: Preinitialized table (input/output argument). */
5142 /* SODITB: Preinitialized table (input/output argument). */
5143 /* DIDITB: Preinitialized table (input/output argument). */
5145 /* OUTPUT ARGUMENTS : */
5146 /* ------------------- */
5147 /* SOSOTB: Table where the terms */
5148 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5149 /* are added with ui and vj positive roots of Legendre polynom */
5150 /* of degree NBPNTU and NBPNTV respectively. */
5151 /* DISOTB: Table where the terms */
5152 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5153 /* are added with ui and vj positive roots of Legendre polynom */
5154 /* of degree NBPNTU and NBPNTV respectively. */
5155 /* SODITB: Table where the terms */
5156 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5157 /* are added with ui and vj positive roots of Legendre polynom */
5158 /* of degree NBPNTU and NBPNTV respectively. */
5159 /* DIDITB: Table where the terms */
5160 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5161 /* are added with ui and vj positive roots of Legendre polynom */
5162 /* of degree NBPNTU and NBPNTV respectively. */
5163 /* FPNTAB: Auxiliary table. */
5164 /* TTABLE: Auxiliary table. */
5165 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5166 /* the returned error code is equal to error code of FONCNP + 100. */
5168 /* COMMONS USED : */
5169 /* ---------------- */
5171 /* REFERENCES CALLED : */
5172 /* --------------------- */
5174 /* DESCRIPTION/NOTES/LIMITATIONS : */
5175 /* ----------------------------------- */
5176 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5177 /* where MA2FXK should be in the following form : */
5178 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,ISOFAV,TCONST,NBPTAB */
5179 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5180 /* with the following input arguments : */
5181 /* - NDIMEN is integer defined as the sum of dimensions of */
5182 /* sub-spaces (i.e. total dimension of the problem). */
5183 /* - UINTFN(2) is a table of 2 reals containing the interval */
5184 /* by u where the function to be approximated is defined */
5185 /* (so it is equal to UIFONC). */
5186 /* - VINTFN(2) is a table of 2 reals containing the interval */
5187 /* by v where the function to be approximated is defined */
5188 /* (so it is equal to VIFONC). */
5189 /* - ISOFAV, is 1 if it is necessary to calculate points with constant u, */
5190 /* is 2 if it is necessary to calculate points with constant v. */
5191 /* Any other value is an error. */
5192 /* - TCONST, real, value of the fixed parameter. Takes values */
5193 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5194 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5195 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5196 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5197 /* 'free' parameter of discretization (v if IISOFAV=1, */
5198 /* u if IISOFAV=2). */
5199 /* - IDERIU, integer, takes values between 0 (position) */
5200 /* and IORDRE(1) (partial derivative of the function by u */
5201 /* of order IORDRE(1) if IORDRE(1) > 0). */
5202 /* - IDERIV, integer, takes values between 0 (position) */
5203 /* and IORDRE(2) (partial derivative of the function by v */
5204 /* of order IORDRE(2) if IORDRE(2) > 0). */
5205 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5206 /* points of the derivative : */
5213 /* and the output arguments aret : */
5214 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5215 /* NBPTAB points calculated in FONCNP. */
5216 /* - IERCOD is, at output the error code of FONCNP. This code */
5217 /* (integer) should be strictly positive if there is a problem. */
5219 /* The input arguments SHOULD NOT be modified under FONCNP.
5222 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5223 /* values of UROOTB and VROOTB are consequently modified. */
5225 /* -->The results of discretisation are ranked in 4 tables */
5226 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5227 /* during the calculation of coefficients of the polynom of approximation. */
5229 /* When NBPNTU is uneven : */
5230 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5231 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5232 /* When NBPNTV is uneven : */
5233 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5234 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5235 /* When NBPNTU and NBPNTV are uneven : */
5236 /* term SOSOTB(0,0) contains F(0,0). */
5239 /* **********************************************************************
5241 /* Name of the routine */
5244 /* --------------------------- Initialization --------------------------
5247 /* Parameter adjustments */
5248 fpntab_dim1 = *ndimen;
5249 fpntab_offset = fpntab_dim1 + 1;
5250 fpntab -= fpntab_offset;
5254 diditb_dim1 = *nbpntu / 2 + 1;
5255 diditb_dim2 = *nbpntv / 2 + 1;
5256 diditb_offset = diditb_dim1 * diditb_dim2;
5257 diditb -= diditb_offset;
5258 soditb_dim1 = *nbpntu / 2;
5259 soditb_dim2 = *nbpntv / 2;
5260 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5261 soditb -= soditb_offset;
5262 disotb_dim1 = *nbpntu / 2;
5263 disotb_dim2 = *nbpntv / 2;
5264 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5265 disotb -= disotb_offset;
5266 sosotb_dim1 = *nbpntu / 2 + 1;
5267 sosotb_dim2 = *nbpntv / 2 + 1;
5268 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5269 sosotb -= sosotb_offset;
5274 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5276 AdvApp2Var_SysBase::mgenmsg_("MMA2DS1", 7L);
5279 if (*isofav < 1 || *isofav > 2) {
5285 /* **********************************************************************
5287 /* --------- Discretization by U on the roots of the polynom of ------ */
5288 /* --------------- Legendre of degree NBPNTU, iso-V by iso-V --------- */
5289 /* **********************************************************************
5293 mma2ds2_(ndimen, &uintfn[1], &vintfn[1], foncnp, nbpntu, nbpntv, &
5294 urootb[1], &vrootb[1], &iuouv, &sosotb[sosotb_offset], &
5295 disotb[disotb_offset], &soditb[soditb_offset], &diditb[
5296 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5298 /* ******************************************************************
5300 /* --------- Discretization by V on the roots of the polynom of ------ */
5301 /* --------------- Legendre of degree NBPNTV, iso-V by iso-V --------- */
5302 /* ******************************************************************
5306 /* --> Inversion of indices of tables */
5308 for (nd = 1; nd <= i__1; ++nd) {
5309 isz1 = *nbpntu / 2 + 1;
5310 isz2 = *nbpntv / 2 + 1;
5311 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5312 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5313 ibid1, &ibid2, iercod);
5317 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5318 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5319 ibid1, &ibid2, iercod);
5325 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5326 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5327 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5331 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5332 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5333 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5340 mma2ds2_(ndimen, &vintfn[1], &uintfn[1], foncnp, nbpntv, nbpntu, &
5341 vrootb[1], &urootb[1], &iuouv, &sosotb[sosotb_offset], &
5342 soditb[soditb_offset], &disotb[disotb_offset], &diditb[
5343 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5344 /* --> Inversion of indices of tables */
5346 for (nd = 1; nd <= i__1; ++nd) {
5347 isz1 = *nbpntv / 2 + 1;
5348 isz2 = *nbpntu / 2 + 1;
5349 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5350 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5351 ibid1, &ibid2, iercod);
5355 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5356 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5357 ibid1, &ibid2, iercod);
5363 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5364 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5365 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5369 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5370 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5371 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5379 /* ------------------------------ The end -------------------------------
5385 AdvApp2Var_SysBase::maermsg_("MMA2DS1", iercod, 7L);
5388 AdvApp2Var_SysBase::mgsomsg_("MMA2DS1", 7L);
5393 //=======================================================================
5394 //function : mma2ds2_
5396 //=======================================================================
5397 int mma2ds2_(integer *ndimen,
5427 static integer c__0 = 0;
5428 /* System generated locals */
5429 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5430 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5431 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5432 fpntab_offset, i__1, i__2, i__3;
5434 /* Local variables */
5435 static integer jdec;
5436 static logical ldbg;
5437 static doublereal alinu, blinu, alinv, blinv, tcons;
5438 static doublereal dbfn1[2], dbfn2[2];
5439 static integer nuroo, nvroo, id, iu, iv;
5440 static doublereal um, up;
5443 /* **********************************************************************
5448 /* Discretization of function F(u,v) on the roots of polynoms of Legendre. */
5452 /* FONCTION&,DISCRETISATION,&POINT */
5454 /* INPUT ARGUMENTS : */
5455 /* ------------------ */
5456 /* NDIMEN: Dimension of the space. */
5457 /* UINTFN: Limits of the interval of definition by u of the function */
5458 /* to be processed: (UINTFN(1),UINTFN(2)). */
5459 /* VINTFN: Limits of the interval of definition by v of the function */
5460 /* to be processed: (VINTFN(1),VINTFN(2)). */
5461 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5462 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5463 /* FONCNP is discretized by u. */
5464 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5465 /* FONCNP is discretized by v. */
5466 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5467 /* of Legendre of degree NBPNTU defined on (-1,1). */
5468 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5469 /* of Legendre of degree NBPNTV defined on (-1,1). */
5470 /* IIUOUV: Shows the type of iso of F(u,v) tom be extracted to improve the */
5471 /* rapidity of calculation (has no influence on the form of result) */
5472 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5473 /* with fixed u (so with NBPNTV values different from v). */
5474 /* = 2, shows that it is necessary to calculate the points of F(u,v) */
5475 /* with fixed v (so with NBPNTV values different from u). */
5476 /* SOSOTB: Preinitialized table (input/output argument). */
5477 /* DISOTB: Preinitialized table (input/output argument). */
5478 /* SODITB: Preinitialized table (input/output argument). */
5479 /* DIDITB: Preinitialized table (input/output argument). */
5481 /* OUTPUT ARGUMENTS : */
5482 /* ------------------- */
5483 /* SOSOTB: Table where the terms */
5484 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5485 /* are added with ui and vj positive roots of Legendre polynom */
5486 /* of degree NBPNTU and NBPNTV respectively. */
5487 /* DISOTB: Table where the terms */
5488 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5489 /* are added with ui and vj positive roots of Legendre polynom */
5490 /* of degree NBPNTU and NBPNTV respectively. */
5491 /* SODITB: Table where the terms */
5492 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5493 /* are added with ui and vj positive roots of Legendre polynom */
5494 /* of degree NBPNTU and NBPNTV respectively. */
5495 /* DIDITB: Table where the terms */
5496 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5497 /* are added with ui and vj positive roots of Legendre polynom */
5498 /* of degree NBPNTU and NBPNTV respectively. */
5499 /* FPNTAB: Auxiliary table. */
5500 /* TTABLE: Auxiliary table. */
5501 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5502 /* the returned error code is equal to error code of FONCNP + 100. */
5504 /* COMMONS USED : */
5505 /* ---------------- */
5507 /* REFERENCES CALLED : */
5508 /* --------------------- */
5510 /* DESCRIPTION/NOTES/LIMITATIONS : */
5511 /* ----------------------------------- */
5512 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5513 /* where MA2FXK should be in the following form : */
5514 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIIUOUV,TCONST,NBPTAB */
5515 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5516 /* with the following input arguments : */
5517 /* - NDIMEN is integer defined as the sum of dimensions of */
5518 /* sub-spaces (i.e. total dimension of the problem). */
5519 /* - UINTFN(2) is a table of 2 reals containing the interval */
5520 /* by u where the function to be approximated is defined */
5521 /* (so it is equal to UIFONC). */
5522 /* - VINTFN(2) is a table of 2 reals containing the interval */
5523 /* by v where the function to be approximated is defined */
5524 /* (so it is equal to VIFONC). */
5525 /* - IIIUOUV, is 1 if it is necessary to calculate points with constant u, */
5526 /* is 2 if it is necessary to calculate points with constant v. */
5527 /* Any other value is an error. */
5528 /* - TCONST, real, value of the fixed parameter. Takes values */
5529 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5530 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5531 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5532 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5533 /* 'free' parameter of discretization (v if IIIUOUV=1, */
5534 /* u if IIIUOUV=2). */
5535 /* - IDERIU, integer, takes values between 0 (position) */
5536 /* and IORDRE(1) (partial derivative of the function by u */
5537 /* of order IORDRE(1) if IORDRE(1) > 0). */
5538 /* - IDERIV, integer, takes values between 0 (position) */
5539 /* and IORDRE(2) (partial derivative of the function by v */
5540 /* of order IORDRE(2) if IORDRE(2) > 0). */
5541 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5542 /* points of the derivative : */
5549 /* and the output arguments aret : */
5550 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5551 /* NBPTAB points calculated in FONCNP. */
5552 /* - IERCOD is, at output the error code of FONCNP. This code */
5553 /* (integer) should be strictly positive if there is a problem. */
5555 /* The input arguments SHOULD NOT be modified under FONCNP.
5558 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5559 /* values of UROOTB and VROOTB are consequently modified. */
5561 /* -->The results of discretisation are ranked in 4 tables */
5562 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5563 /* during the calculation of coefficients of the polynom of approximation. */
5565 /* When NBPNTU is uneven : */
5566 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5567 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5568 /* When NBPNTV is uneven : */
5569 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5570 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5571 /* When NBPNTU and NBPNTV are uneven : */
5572 /* term SOSOTB(0,0) contains F(0,0). */
5574 /* ATTENTION: These 4 tables are filled by varying the */
5575 /* 1st index first. So, the discretizations */
5576 /* of F(...,t) (for IIUOUV = 2) or of F(t,...) (IIUOUV = 1) */
5577 /* are stored in SOSOTB(...,t), SODITB(...,t), etc... */
5578 /* (this allows to gain important time). */
5579 /* It is required that the caller, in case of IIUOUV=1, */
5580 /* invert the roles of u and v, of SODITB and DISOTB BEFORE the */
5583 /* **********************************************************************
5586 /* Name of the routine */
5588 /* --> Indices of loops. */
5590 /* --------------------------- Initialization --------------------------
5593 /* Parameter adjustments */
5597 fpntab_dim1 = *ndimen;
5598 fpntab_offset = fpntab_dim1 + 1;
5599 fpntab -= fpntab_offset;
5601 diditb_dim1 = *nbpntu / 2 + 1;
5602 diditb_dim2 = *nbpntv / 2 + 1;
5603 diditb_offset = diditb_dim1 * diditb_dim2;
5604 diditb -= diditb_offset;
5605 soditb_dim1 = *nbpntu / 2;
5606 soditb_dim2 = *nbpntv / 2;
5607 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5608 soditb -= soditb_offset;
5609 disotb_dim1 = *nbpntu / 2;
5610 disotb_dim2 = *nbpntv / 2;
5611 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5612 disotb -= disotb_offset;
5613 sosotb_dim1 = *nbpntu / 2 + 1;
5614 sosotb_dim2 = *nbpntv / 2 + 1;
5615 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5616 sosotb -= sosotb_offset;
5620 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5622 AdvApp2Var_SysBase::mgenmsg_("MMA2DS2", 7L);
5626 alinu = (uintfn[2] - uintfn[1]) / 2.;
5627 blinu = (uintfn[2] + uintfn[1]) / 2.;
5628 alinv = (vintfn[2] - vintfn[1]) / 2.;
5629 blinv = (vintfn[2] + vintfn[1]) / 2.;
5632 dbfn1[0] = vintfn[1];
5633 dbfn1[1] = vintfn[2];
5634 dbfn2[0] = uintfn[1];
5635 dbfn2[1] = uintfn[2];
5637 dbfn1[0] = uintfn[1];
5638 dbfn1[1] = uintfn[2];
5639 dbfn2[0] = vintfn[1];
5640 dbfn2[1] = vintfn[2];
5643 /* **********************************************************************
5645 /* -------- Discretization by U on the roots of Legendre polynom -------- */
5646 /* ---------------- of degree NBPNTU, with Vj fixed -------------------- */
5647 /* **********************************************************************
5650 nuroo = *nbpntu / 2;
5651 nvroo = *nbpntv / 2;
5652 jdec = (*nbpntu + 1) / 2;
5654 /* ----------- Loading of parameters of discretization by U ------------- */
5657 for (iu = 1; iu <= i__1; ++iu) {
5658 ttable[iu] = blinu + alinu * urootb[iu];
5662 /* -------------- For Vj fixed, negative root of Legendre ------------- */
5665 for (iv = 1; iv <= i__1; ++iv) {
5666 tcons = blinv + alinv * vrootb[iv];
5667 (*foncnp)(ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5668 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5673 for (id = 1; id <= i__2; ++id) {
5675 for (iu = 1; iu <= i__3; ++iu) {
5676 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5677 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5678 sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1]
5679 = sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) *
5680 sosotb_dim1] + up + um;
5681 disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) * disotb_dim1]
5682 = disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) *
5683 disotb_dim1] + up - um;
5684 soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) * soditb_dim1]
5685 = soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) *
5686 soditb_dim1] - up - um;
5687 diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1]
5688 = diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) *
5689 diditb_dim1] - up + um;
5692 if (*nbpntu % 2 != 0) {
5693 up = fpntab[id + jdec * fpntab_dim1];
5694 sosotb[(nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1] +=
5696 diditb[(nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1] -=
5704 /* --------- For Vj = 0 (uneven NBPNTV), discretization by U ----------- */
5706 if (*nbpntv % 2 != 0) {
5708 (*foncnp)(ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5709 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5714 for (id = 1; id <= i__1; ++id) {
5716 for (iu = 1; iu <= i__2; ++iu) {
5717 up = fpntab[id + (jdec + iu) * fpntab_dim1];
5718 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5719 sosotb[iu + id * sosotb_dim2 * sosotb_dim1] = sosotb[iu + id *
5720 sosotb_dim2 * sosotb_dim1] + up + um;
5721 diditb[iu + id * diditb_dim2 * diditb_dim1] = diditb[iu + id *
5722 diditb_dim2 * diditb_dim1] + up - um;
5725 if (*nbpntu % 2 != 0) {
5726 up = fpntab[id + jdec * fpntab_dim1];
5727 sosotb[id * sosotb_dim2 * sosotb_dim1] += up;
5733 /* -------------- For Vj fixed, positive root of Legendre ------------- */
5736 for (iv = 1; iv <= i__1; ++iv) {
5737 tcons = alinv * vrootb[(*nbpntv + 1) / 2 + iv] + blinv;
5738 (*foncnp)(ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5739 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5744 for (id = 1; id <= i__2; ++id) {
5746 for (iu = 1; iu <= i__3; ++iu) {
5747 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5748 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5749 sosotb[iu + (iv + id * sosotb_dim2) * sosotb_dim1] = sosotb[
5750 iu + (iv + id * sosotb_dim2) * sosotb_dim1] + up + um;
5751 disotb[iu + (iv + id * disotb_dim2) * disotb_dim1] = disotb[
5752 iu + (iv + id * disotb_dim2) * disotb_dim1] + up - um;
5753 soditb[iu + (iv + id * soditb_dim2) * soditb_dim1] = soditb[
5754 iu + (iv + id * soditb_dim2) * soditb_dim1] + up + um;
5755 diditb[iu + (iv + id * diditb_dim2) * diditb_dim1] = diditb[
5756 iu + (iv + id * diditb_dim2) * diditb_dim1] + up - um;
5759 if (*nbpntu % 2 != 0) {
5760 up = fpntab[id + jdec * fpntab_dim1];
5761 sosotb[(iv + id * sosotb_dim2) * sosotb_dim1] += up;
5762 diditb[(iv + id * diditb_dim2) * diditb_dim1] += up;
5769 /* ------------------------------ The end -------------------------------
5775 AdvApp2Var_SysBase::maermsg_("MMA2DS2", iercod, 7L);
5778 AdvApp2Var_SysBase::mgsomsg_("MMA2DS2", 7L);
5783 //=======================================================================
5784 //function : mma2er1_
5786 //=======================================================================
5787 int mma2er1_(integer *ndjacu,
5803 /* System generated locals */
5804 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
5807 /* Local variables */
5808 static logical ldbg;
5809 static integer minu, minv;
5810 static doublereal vaux[2];
5811 static integer ii, nd, jj;
5812 static doublereal bid0, bid1;
5815 /* **********************************************************************
5820 /* Calculate max approximation error done when */
5821 /* the coefficients of PATJAC such that the degree by U varies between */
5822 /* MINDGU and MAXDGU and the degree by V varies between MINDGV and MAXDGV are removed. */
5826 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5828 /* INPUT ARGUMENTS : */
5829 /* ------------------ */
5830 /* NDJACU: Dimension by U of table PATJAC. */
5831 /* NDJACV: Dimension by V of table PATJAC. */
5832 /* NDIMEN: Dimension of the space. */
5833 /* MINDGU: Lower limit of index by U of coeff. of PATJAC to be taken into account. */
5834 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5835 /* MINDGV: Lower limit of index by V of coeff. of PATJAC to be taken into account. */
5836 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5837 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5838 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5839 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5840 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5841 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5842 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5843 /* PATJAC: Table of coeff. of square of approximation with */
5844 /* constraints of order IORDRU by U and IORDRV by V. */
5845 /* VECERR: Auxiliary vector. */
5846 /* ERREUR: MAX Error commited during removal of ALREADY CALCULATED coeff of PATJAC */
5848 /* OUTPUT ARGUMENTS : */
5849 /* ------------------- */
5850 /* ERREUR: MAX Error commited during removal of coeff of PATJAC */
5851 /* of indices from MINDGU to MAXDGU by U and from MINDGV to MAXDGV by V */
5852 /* THEN the already calculated error. */
5854 /* COMMONS USED : */
5855 /* ---------------- */
5857 /* REFERENCES CALLED : */
5858 /* --------------------- */
5860 /* DESCRIPTION/NOTES/LIMITATIONS : */
5861 /* ----------------------------------- */
5862 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5863 /* approximation of F(U,V). The indices i and j show the degree */
5864 /* by U and by V of base polynoms. These polynoms have the form: */
5866 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5868 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5869 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5871 /* The contribution to the error of term Cij when it is */
5872 /* removed from PATJAC is increased by: */
5874 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5876 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5878 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5882 /* ***********************************************************************
5884 /* Name of the routine */
5887 /* ----------------------------- Initialisations ------------------------
5890 /* Parameter adjustments */
5892 patjac_dim1 = *ndjacu + 1;
5893 patjac_dim2 = *ndjacv + 1;
5894 patjac_offset = patjac_dim1 * patjac_dim2;
5895 patjac -= patjac_offset;
5898 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5900 AdvApp2Var_SysBase::mgenmsg_("MMA2ER1", 7L);
5903 minu = (*iordru + 1) << 1;
5904 minv = (*iordrv + 1) << 1;
5906 /* ------------------- Calculate the increment of the max error --------------- */
5907 /* ----- during the removal of the coeffs of indices from MINDGU to MAXDGU ---- */
5908 /* ---------------- by U and indices from MINDGV to MAXDGV by V --------------- */
5911 for (nd = 1; nd <= i__1; ++nd) {
5914 for (jj = *mindgv; jj <= i__2; ++jj) {
5917 for (ii = *mindgu; ii <= i__3; ++ii) {
5918 bid0 += (d__1 = patjac[ii + (jj + nd * patjac_dim2) *
5919 patjac_dim1], abs(d__1)) * xmaxju[ii - minu];
5922 bid1 = bid0 * xmaxjv[jj - minv] + bid1;
5930 /* ----------------------- Calculate the max error ----------------------*/
5932 bid1 = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
5936 *erreur = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
5938 /* ------------------------- The end ------------------------------------
5942 AdvApp2Var_SysBase::mgsomsg_("MMA2ER1", 7L);
5947 //=======================================================================
5948 //function : mma2er2_
5950 //=======================================================================
5951 int mma2er2_(integer *ndjacu,
5963 doublereal *epmscut,
5970 /* System generated locals */
5971 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2;
5974 /* Local variables */
5975 static logical ldbg;
5976 static doublereal vaux[2];
5977 static integer i2rdu, i2rdv;
5978 static doublereal errnu, errnv;
5979 static integer ii, nd, jj, nu, nv;
5980 static doublereal bid0, bid1;
5983 /* **********************************************************************
5988 /* Remove coefficients of PATJAC to obtain the minimum degree */
5989 /* by U and V checking the imposed tolerance. */
5993 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5995 /* INPUT ARGUMENTS : */
5996 /* ------------------ */
5997 /* NDJACU: Degree by U of table PATJAC. */
5998 /* NDJACV: Degree by V of table PATJAC. */
5999 /* NDIMEN: Dimension of the space. */
6000 /* MINDGU: Limit of index by U of coeff. of PATJAC to be PRESERVED (should be >=0). */
6001 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
6002 /* MINDGV: Limit of index by V of coeff. of PATJAC to be PRESERVED (should be >=0). */
6003 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
6004 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
6005 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
6006 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
6007 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
6008 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
6009 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
6010 /* PATJAC: Table of coeff. of square of approximation with */
6011 /* constraints of order IORDRU by U and IORDRV by V. */
6012 /* EPMSCUT: Tolerance of approximation. */
6013 /* VECERR: Auxiliary vector. */
6014 /* ERREUR: MAX Error commited ALREADY CALCULATED */
6016 /* OUTPUT ARGUMENTS : */
6017 /* ------------------- */
6018 /* ERREUR: MAX Error commited by preserving only coeff of PATJAC */
6019 /* of indices from 0 to NEWDGU by U and from 0 to NEWDGV by V */
6020 /* PLUS the already calculated error. */
6021 /* NEWDGU: Min. Degree by U such as the square of approximation */
6022 /* could check the tolerance. There is always NEWDGU >= MINDGU >= 0. */
6023 /* NEWDGV: Min. Degree by V such as the square of approximation */
6024 /* could check the tolerance. There is always NEWDGV >= MINDGV >= 0. */
6027 /* COMMONS USED : */
6028 /* ---------------- */
6030 /* REFERENCES CALLED : */
6031 /* --------------------- */
6033 /* DESCRIPTION/NOTES/LIMITATIONS : */
6034 /* ----------------------------------- */
6035 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
6036 /* approximation of F(U,V). The indices i and j show the degree */
6037 /* by U and by V of base polynoms. These polynoms have the form: */
6039 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
6041 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
6042 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
6044 /* The contribution to the error of term Cij when it is */
6045 /* removed from PATJAC is increased by: */
6047 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
6049 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
6051 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
6055 /* **********************************************************************
6057 /* Name of the routine */
6060 /* ----------------------------- Initialisations ------------------------
6063 /* Parameter adjustments */
6065 patjac_dim1 = *ndjacu + 1;
6066 patjac_dim2 = *ndjacv + 1;
6067 patjac_offset = patjac_dim1 * patjac_dim2;
6068 patjac -= patjac_offset;
6071 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
6073 AdvApp2Var_SysBase::mgenmsg_("MMA2ER2", 7L);
6076 i2rdu = (*iordru + 1) << 1;
6077 i2rdv = (*iordrv + 1) << 1;
6081 /* **********************************************************************
6083 /* -------------------- Cutting of oefficients ------------------------
6085 /* **********************************************************************
6090 /* ------------------- Calculate the increment of max error --------------- */
6091 /* ----- during the removal of coeff. of indices from MINDGU to MAXDGU ------ */
6092 /* ---------------- by U, the degree by V is fixed to NV -----------------
6097 bid0 = xmaxjv[nv - i2rdv];
6099 for (nd = 1; nd <= i__1; ++nd) {
6102 for (ii = i2rdu; ii <= i__2; ++ii) {
6103 bid1 += (d__1 = patjac[ii + (nv + nd * patjac_dim2) *
6104 patjac_dim1], abs(d__1)) * xmaxju[ii - i2rdu] * bid0;
6111 vecerr[1] = *epmscut * 2;
6113 errnv = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6115 /* ------------------- Calculate the increment of max error --------------- */
6116 /* ----- during the removal of coeff. of indices from MINDGV to MAXDGV ------ */
6117 /* ---------------- by V, the degree by U is fixed to NU -----------------
6122 bid0 = xmaxju[nu - i2rdu];
6124 for (nd = 1; nd <= i__1; ++nd) {
6127 for (jj = i2rdv; jj <= i__2; ++jj) {
6128 bid1 += (d__1 = patjac[nu + (jj + nd * patjac_dim2) *
6129 patjac_dim1], abs(d__1)) * xmaxjv[jj - i2rdv] * bid0;
6136 vecerr[1] = *epmscut * 2;
6138 errnu = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6140 /* ----------------------- Calculate the max error ----------------------
6146 errnu = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6148 errnv = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6150 if (errnu > errnv) {
6151 if (errnv < *epmscut) {
6158 if (errnu < *epmscut) {
6168 /* -------------------------- Return the degrees -------------------
6172 *newdgu = max(nu,1);
6173 *newdgv = max(nv,1);
6175 /* ----------------------------------- The end --------------------------
6179 AdvApp2Var_SysBase::mgsomsg_("MMA2ER2", 7L);
6184 //=======================================================================
6185 //function : mma2fnc_
6187 //=======================================================================
6188 int AdvApp2Var_ApproxF2var::mma2fnc_(integer *ndimen,
6228 static integer c__8 = 8;
6230 /* System generated locals */
6231 integer courbe_dim1, courbe_dim2, courbe_offset, somtab_dim1, somtab_dim2,
6232 somtab_offset, diftab_dim1, diftab_dim2, diftab_offset,
6233 contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
6234 contr2_offset, errmax_dim1, errmax_offset, errmoy_dim1,
6235 errmoy_offset, i__1;
6238 /* Local variables */
6239 static integer ideb;
6240 static doublereal tmil;
6241 static integer ideb1, ibid1, ibid2, ncfja, ndgre, ilong,
6243 static doublereal wrkar[1];
6244 static integer nupil;
6245 static long int iofwr;
6246 static doublereal uvpav[4] /* was [2][2] */;
6247 static integer nd, ii;
6250 static doublereal uv11[4] /* was [2][2] */;
6251 static integer ncb1;
6252 static doublereal eps3;
6253 static integer isz1, isz2, isz3, isz4, isz5;
6254 static long int ipt1, ipt2, ipt3, ipt4, ipt5,iptt, jptt;
6256 /* **********************************************************************
6261 /* Approximation of a limit of non polynomial function F(u,v) */
6262 /* (in the space of dimension NDIMEN) by SEVERAL */
6263 /* polynomial curves, by the method of least squares. The parameter of the function is preserved. */
6267 /* TOUS, AB_SPECIFI :: FONCTION&,EXTREMITE&, APPROXIMATION, &COURBE. */
6269 /* INPUT ARGUMENTS : */
6270 /* ----------------- */
6271 /* NDIMEN: Total Dimension of the space (sum of dimensions */
6272 /* of sub-spaces) */
6273 /* NBSESP: Number of "independent" sub-spaces. */
6274 /* NDIMSE: Table of dimensions of sub-spaces. */
6275 /* UVFONC: Limits of the interval (a,b)x(c,d) of definition of the */
6276 /* function to be approached by U (UVFONC(*,1) contains (a,b)) */
6277 /* and by V (UVFONC(*,2) contains (c,d)). */
6278 /* FONCNP: External function of position on the non polynomial function to be approached. */
6279 /* TCONST: Value of isoparameter of F(u,v) to be discretized. */
6280 /* ISOFAV: Type of chosen iso, = 1, shose that discretization is with u */
6281 /* fixed; = 2, shows that v is fixed. */
6282 /* NBROOT: Nb of points of discretisation of the iso, extremities not included. */
6283 /* ROOTLG: Table of roots of the polynom of Legendre defined on */
6284 /* (-1,1), of degree NBROOT. */
6285 /* IORDRE: Order of constraint at the extremities of the limit */
6286 /* -1 = no constraints, */
6287 /* 0 = constraints of passage to limits (i.e. C0), */
6288 /* 1 = C0 + constraints of 1st derivatives (i.e. C1), */
6289 /* 2 = C1 + constraints of 2nd derivatives (i.e. C2). */
6290 /* IDERIV: Order of derivative of the limit. */
6291 /* NDGJAC: Degree of serial development to be used for calculation in */
6292 /* the Jacobi base. */
6293 /* NBCRMX: Max Nb of curves to be created. */
6294 /* NCFLIM: Max Nb of coeff of the polynomial curve */
6295 /* of approximation (should be above or equal to */
6296 /* 2*IORDRE+2 and below or equal to 50). */
6297 /* EPSAPR: Table of required errors of approximation */
6298 /* sub-space by sub-space. */
6300 /* OUTPUT ARGUMENTS : */
6301 /* ------------------- */
6302 /* NCOEFF: Number of significative coeff of calculated curves. */
6303 /* COURBE: Table of coeff. of calculated polynomial curves. */
6304 /* Should be dimensioned in (NCFLIM,NDIMEN,NBCRMX). */
6305 /* These curves are ALWAYS parametrized in (-1,1). */
6306 /* NBCRBE: Nb of calculated curves. */
6307 /* SOMTAB: For F defined on (-1,1) (otherwise rescale the */
6308 /* parameters), this is the table of sums F(u,vj) + F(u,-vj)
6310 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6311 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6312 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6313 /* v of order IDERIV are taken). */
6314 /* DIFTAB: For F defined on (-1,1) (otherwise rescale the */
6315 /* parameters), this is the table of sums F(u,vj) - F(u,-vj)
6317 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6318 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6319 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6320 /* v of order IDERIV are taken). */
6321 /* CONTR1: Contains the coordinates of the left extremity of the iso */
6322 /* and of its derivatives till order IORDRE */
6323 /* CONTR2: Contains the coordinates of the right extremity of the iso */
6324 /* and of its derivatives till order IORDRE */
6325 /* TABDEC: Table of NBCRBE+1 parameters of cut of UVFONC(1:2,1)
6327 /* if ISOFAV=2, or of UVFONC(1:2,2) if ISOFAV=1. */
6328 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6329 /* committed in the approximation of FONCNP by NBCRBE curves. */
6330 /* ERRMOY: Table of AVERAGE errors (sub-space by sub-space) */
6331 /* committed in the approximation of FONCNP by NBCRBE curves.
6332 /* IERCOD: Error code: */
6333 /* -1 = ERRMAX > EPSAPR for at least one sub-space. */
6334 /* (the resulting curves of at least mathematic degree NCFLIM-1 */
6335 /* are calculated). */
6336 /* 0 = Everything is ok. */
6337 /* 1 = Pb of incoherence of inputs. */
6338 /* 10 = Pb of calculation of the interpolation of constraints. */
6339 /* 13 = Pb in the dynamic allocation. */
6340 /* 33 = Pb in the data recuperation from block data */
6341 /* of coeff. of integration by GAUSS method. */
6342 /* >100 Pb in the evaluation of FONCNP, the returned error code */
6343 /* is equal to the error code of FONCNP + 100. */
6345 /* COMMONS USED : */
6346 /* ---------------- */
6348 /* REFERENCES CALLED : */
6349 /* ----------------------- */
6351 /* DESCRIPTION/NOTES/LIMITATIONS : */
6352 /* ----------------------------------- */
6353 /* --> The approximation part is done in the space of dimension */
6354 /* NDIMEN (the sum of dimensions of sub-spaces). For example : */
6355 /* If NBSESP=2 and NDIMSE(1)=3, NDIMSE(2)=2, there is smoothing with */
6356 /* NDIMEN=5. The result (in COURBE(NDIMEN,NCOEFF,i) ), will be */
6357 /* composed of the result of smoothing of 3D function in */
6358 /* COURBE(1:3,1:NCOEFF,i) and of smoothing of 2D function in */
6359 /* COURBE(4:5,1:NCOEFF,i). */
6361 /* --> Routine FONCNP should be declared EXTERNAL in the program */
6362 /* calling MMA2FNC. */
6364 /* --> Function FONCNP, declared externally, should be declared */
6365 /* IMPERATIVELY in form : */
6366 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIUOUV,TCONST,NBPTAB */
6367 /* ,TTABLE,IDERIU,IDERIV,IERCOD) */
6368 /* where the input arguments are : */
6369 /* - NDIMEN is integer defined as the sum of dimensions of */
6370 /* sub-spaces (i.e. total dimension of the problem). */
6371 /* - UINTFN(2) is a table of 2 reals containing the interval */
6372 /* by u where the function to be approximated is defined */
6373 /* (so it is equal to UIFONC). */
6374 /* - VINTFN(2) is a table of 2 reals containing the interval */
6375 /* by v where the function to be approximated is defined */
6376 /* (so it is equal to VIFONC). */
6377 /* - IIUOUV, shows that the points to be calculated have a constant U */
6378 /* (IIUOUV=1) or a constant V (IIUOUV=2). */
6379 /* - TCONST, real, value of the fixed discretisation parameter. Takes values */
6380 /* in (UINTFN(1),UINTFN(2)) if IIUOUV=1, */
6381 /* or in (VINTFN(1),VINTFN(2)) if IIUOUV=2. */
6382 /* - NBPTAB, the nb of point of discretisation following the free variable */
6383 /* : V if IIUOUV=1 or U if IIUOUV = 2. */
6384 /* - TTABLE, Table of NBPTAB parametres of discretisation. . */
6385 /* - IDERIU, integer, takes values between 0 (position) */
6386 /* and IORDREU (partial derivative of the function by u */
6387 /* of order IORDREU if IORDREU > 0). */
6388 /* - IDERIV, integer, takes values between 0 (position) */
6389 /* and IORDREV (partial derivative of the function by v */
6390 /* of order IORDREV if IORDREV > 0). */
6391 /* and the output arguments are : */
6392 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
6393 /* NBPTAB points calculated in FONCNP. */
6394 /* - IERCOD is, at output the error code of FONCNP. This code */
6395 /* (integer) should be strictly positive if there is a problem. */
6397 /* The input arguments SHOULD NOT BE modified under FONCNP.
6400 /* --> If IERCOD=-1, the required precision can't be reached (ERRMAX */
6401 /* is above EPSAPR on at least one sub-space), but
6403 /* one gives the best possible result for NCFLIM and EPSAPR */
6404 /* chosen by the user. In this case (and for IERCOD=0), there is a solution. */
6407 /* **********************************************************************
6409 /* Name of the routine */
6411 /* Parameter adjustments */
6416 errmoy_dim1 = *nbsesp;
6417 errmoy_offset = errmoy_dim1 + 1;
6418 errmoy -= errmoy_offset;
6419 errmax_dim1 = *nbsesp;
6420 errmax_offset = errmax_dim1 + 1;
6421 errmax -= errmax_offset;
6422 contr2_dim1 = *ndimen;
6423 contr2_dim2 = *iordre + 2;
6424 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
6425 contr2 -= contr2_offset;
6426 contr1_dim1 = *ndimen;
6427 contr1_dim2 = *iordre + 2;
6428 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
6429 contr1 -= contr1_offset;
6430 diftab_dim1 = *nbroot / 2 + 1;
6431 diftab_dim2 = *ndimen;
6432 diftab_offset = diftab_dim1 * (diftab_dim2 + 1);
6433 diftab -= diftab_offset;
6434 somtab_dim1 = *nbroot / 2 + 1;
6435 somtab_dim2 = *ndimen;
6436 somtab_offset = somtab_dim1 * (somtab_dim2 + 1);
6437 somtab -= somtab_offset;
6439 courbe_dim1 = *ncflim;
6440 courbe_dim2 = *ndimen;
6441 courbe_offset = courbe_dim1 * (courbe_dim2 + 1) + 1;
6442 courbe -= courbe_offset;
6445 ibb = AdvApp2Var_SysBase::mnfndeb_();
6447 AdvApp2Var_SysBase::mgenmsg_("MMA2FNC", 7L);
6452 /* ---------------- Set to zero the coefficients of CURVE --------------
6455 ilong = *ndimen * *ncflim * *nbcrmx;
6456 AdvApp2Var_SysBase::mvriraz_(&ilong, (char *)&courbe[courbe_offset]);
6458 /* **********************************************************************
6460 /* -------------------------- Checking of entries ------------------
6462 /* **********************************************************************
6465 AdvApp2Var_MathBase::mmveps3_(&eps3);
6466 if ((d__1 = uvfonc[4] - uvfonc[3], abs(d__1)) < eps3) {
6469 if ((d__1 = uvfonc[6] - uvfonc[5], abs(d__1)) < eps3) {
6478 /* ********************************************************************** */
6479 /* ------------- Preparation of parameters of discretisation ----------- */
6480 /* **********************************************************************
6483 /* -- Allocation of a table of parameters and points of discretisation -- */
6484 /* --> For the parameters of discretisation. */
6486 /* --> For the points of discretisation in MMA1FDI and MMA1CDI and
6488 /* the auxiliary curve for MMAPCMP */
6489 ibid1 = *ndimen * (*nbroot + 2);
6490 ibid2 = ((*iordre + 1) << 1) * *nbroot;
6491 isz2 = max(ibid1,ibid2);
6492 ibid1 = (((*ncflim - 1) / 2 + 1) << 1) * *ndimen;
6493 isz2 = max(ibid1,isz2);
6494 /* --> To return the polynoms of hermit. */
6495 isz3 = ((*iordre + 1) << 2) * (*iordre + 1);
6496 /* --> For the Gauss coeff. of integration. */
6497 isz4 = (*nbroot / 2 + 1) * (*ndgjac + 1 - ((*iordre + 1) << 1));
6498 /* --> For the coeff of the curve in the base of Jacobi */
6499 isz5 = (*ndgjac + 1) * *ndimen;
6501 ndwrk = isz1 + isz2 + isz3 + isz4 + isz5;
6502 AdvApp2Var_SysBase::mcrrqst_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6505 /* --> For the parameters of discretisation (NBROOT+2 extremities). */
6507 /* --> For the points of discretisation FPNTAB(NDIMEN,NBROOT+2), */
6508 /* FPNTAB(NBROOT,2*(IORDRE+1)) and for WRKAR of MMAPCMP. */
6510 /* --> For the polynoms of Hermit */
6512 /* --> For the Gauss coeff of integration. */
6514 /* --> For the curve in Jacobi. */
6517 /* ------------------ Initialisation of management of cuts ---------
6521 uvpav[0] = uvfonc[3];
6522 uvpav[1] = uvfonc[4];
6523 tabdec[0] = uvfonc[5];
6524 tabdec[1] = uvfonc[6];
6525 } else if (*isofav == 2) {
6526 tabdec[0] = uvfonc[3];
6527 tabdec[1] = uvfonc[4];
6528 uvpav[2] = uvfonc[5];
6529 uvpav[3] = uvfonc[6];
6537 /* **********************************************************************
6539 /* APPROXIMATION WITH CUTS */
6540 /* **********************************************************************
6544 /* --> When the top is reached, this is the end ! */
6545 if (nupil - *nbcrbe == 0) {
6550 uvpav[2] = tabdec[*nbcrbe];
6551 uvpav[3] = tabdec[*nbcrbe + 1];
6552 } else if (*isofav == 2) {
6553 uvpav[0] = tabdec[*nbcrbe];
6554 uvpav[1] = tabdec[*nbcrbe + 1];
6559 /* -------------------- Normalization of parameters -------------------- */
6561 mma1nop_(nbroot, &rootlg[1], uvpav, isofav, &wrkar[ipt1], &ier);
6566 /* -------------------- Discretisation of FONCNP ------------------------ */
6568 mma1fdi_(ndimen, uvpav, foncnp, isofav, tconst, nbroot, &wrkar[ipt1],
6569 iordre, ideriv, &wrkar[ipt2], &somtab[(ncb1 * somtab_dim2 + 1) *
6570 somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6571 contr1[(ncb1 * contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1
6572 * contr2_dim2 + 1) * contr2_dim1 + 1], iercod);
6577 /* -----------Cut the discretisation of constraints ------------*/
6580 mma1cdi_(ndimen, nbroot, &rootlg[1], iordre, &contr1[(ncb1 *
6581 contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1 *
6582 contr2_dim2 + 1) * contr2_dim1 + 1], &somtab[(ncb1 *
6583 somtab_dim2 + 1) * somtab_dim1], &diftab[(ncb1 * diftab_dim2
6584 + 1) * diftab_dim1], &wrkar[ipt2], &wrkar[ipt3], &ier);
6590 /* **********************************************************************
6592 /* -------------------- Calculate the curve of approximation -------------
6594 /* **********************************************************************
6597 mma1jak_(ndimen, nbroot, iordre, ndgjac, &somtab[(ncb1 * somtab_dim2 + 1)
6598 * somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6599 wrkar[ipt4], &wrkar[ipt5], &ier);
6604 /* **********************************************************************
6606 /* ---------------- Add polynom of interpolation -------------------
6608 /* **********************************************************************
6612 mma1cnt_(ndimen, iordre, &contr1[(ncb1 * contr1_dim2 + 1) *
6613 contr1_dim1 + 1], &contr2[(ncb1 * contr2_dim2 + 1) *
6614 contr2_dim1 + 1], &wrkar[ipt3], ndgjac, &wrkar[ipt5]);
6617 /* **********************************************************************
6619 /* --------------- Calculate Max and Average error ----------------------
6621 /* **********************************************************************
6624 mma1fer_(ndimen, nbsesp, &ndimse[1], iordre, ndgjac, &wrkar[ipt5], ncflim,
6625 &epsapr[1], &wrkar[ipt2], &errmax[ncb1 * errmax_dim1 + 1], &
6626 errmoy[ncb1 * errmoy_dim1 + 1], &ncoeff[ncb1], &ier);
6631 if (ier == 0 || (ier == -1 && nupil == *nbcrmx)) {
6633 /* ******************************************************************
6635 /* ----------------------- Compression du resultat ------------------
6637 /* ******************************************************************
6643 ncfja = *ndgjac + 1;
6644 /* -> Compression of result in WRKAR(IPT2) */
6647 AdvApp2Var_MathBase::mmapcmp_(ndimen,
6648 &ncfja, &ncoeff[ncb1], &wrkar[ipt5], &wrkar[ipt2]);
6650 AdvApp2Var_MathBase::mmapcmp_((integer*)ndimen,
6656 ilong = *ndimen * *ncflim;
6657 AdvApp2Var_SysBase::mvriraz_(&ilong, (char*)&wrkar[ipt5]);
6658 /* -> Passage to canonic base (-1,1) (result in WRKAR(IPT5)).
6660 ndgre = ncoeff[ncb1] - 1;
6662 for (nd = 1; nd <= i__1; ++nd) {
6663 iptt = ipt2 + ((nd - 1) << 1) * (ndgre / 2 + 1);
6664 jptt = ipt5 + (nd - 1) * ncoeff[ncb1];
6665 AdvApp2Var_MathBase::mmjacan_(iordre, &ndgre, &wrkar[iptt], &wrkar[jptt]);
6669 /* -> Store the calculated curve */
6671 AdvApp2Var_MathBase::mmfmca8_(&ncoeff[ncb1], ndimen, &ibid1, ncflim, ndimen, &ibid1, &
6672 wrkar[ipt5], &courbe[(ncb1 * courbe_dim2 + 1) * courbe_dim1 +
6675 /* -> Before normalization of constraints on (-1,1), recalculate */
6676 /* the true constraints. */
6678 for (ii = 0; ii <= i__1; ++ii) {
6679 mma1noc_(uv11, ndimen, &ii, &contr1[(ii + 1 + ncb1 * contr1_dim2)
6680 * contr1_dim1 + 1], uvpav, isofav, ideriv, &contr1[(ii +
6681 1 + ncb1 * contr1_dim2) * contr1_dim1 + 1]);
6682 mma1noc_(uv11, ndimen, &ii, &contr2[(ii + 1 + ncb1 * contr2_dim2)
6683 * contr2_dim1 + 1], uvpav, isofav, ideriv, &contr2[(ii +
6684 1 + ncb1 * contr2_dim2) * contr2_dim1 + 1]);
6688 ibid1 = (*nbroot / 2 + 1) * *ndimen;
6689 mma1noc_(uv11, &ibid1, &ii, &somtab[(ncb1 * somtab_dim2 + 1) *
6690 somtab_dim1], uvpav, isofav, ideriv, &somtab[(ncb1 *
6691 somtab_dim2 + 1) * somtab_dim1]);
6692 mma1noc_(uv11, &ibid1, &ii, &diftab[(ncb1 * diftab_dim2 + 1) *
6693 diftab_dim1], uvpav, isofav, ideriv, &diftab[(ncb1 *
6694 diftab_dim2 + 1) * diftab_dim1]);
6697 for (nd = 1; nd <= i__1; ++nd) {
6698 mma1noc_(uv11, &ncoeff[ncb1], &ii, &courbe[(nd + ncb1 *
6699 courbe_dim2) * courbe_dim1 + 1], uvpav, isofav, ideriv, &
6700 courbe[(nd + ncb1 * courbe_dim2) * courbe_dim1 + 1]);
6704 /* -> Update the nb of already created curves */
6707 /* -> ...otherwise try to cut the current interval in 2... */
6709 tmil = (tabdec[*nbcrbe + 1] + tabdec[*nbcrbe]) / 2.;
6712 ilong = (nupil - *nbcrbe) << 3;
6713 AdvApp2Var_SysBase::mcrfill_(&ilong, (char *)&tabdec[ideb],(char *)&tabdec[ideb1]);
6714 tabdec[ideb] = tmil;
6718 /* ---------- Make approximation of the rest -----------
6723 /* --------------------- Return code of error -----------------
6725 /* --> Pb with dynamic allocation */
6729 /* --> Inputs incoherent. */
6734 /* -------------------------- Dynamic desallocation -------------------
6739 AdvApp2Var_SysBase::mcrdelt_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6746 /* ------------------------------ The end -------------------------------
6751 AdvApp2Var_SysBase::maermsg_("MMA2FNC", iercod, 7L);
6754 AdvApp2Var_SysBase::mgsomsg_("MMA2FNC", 7L);
6759 //=======================================================================
6760 //function : mma2fx6_
6762 //=======================================================================
6763 int AdvApp2Var_ApproxF2var::mma2fx6_(integer *ncfmxu,
6780 /* System generated locals */
6781 integer epsfro_dim1, epsfro_offset, patcan_dim1, patcan_dim2, patcan_dim3,
6782 patcan_dim4, patcan_offset, errmax_dim1, errmax_dim2,
6783 errmax_offset, ncoefu_dim1, ncoefu_offset, ncoefv_dim1,
6784 ncoefv_offset, i__1, i__2, i__3, i__4, i__5;
6785 doublereal d__1, d__2;
6787 /* Local variables */
6788 static integer idim, ncfu, ncfv, id, ii, nd, jj, ku, kv, ns, ibb;
6789 static doublereal bid;
6790 static doublereal tol;
6792 /* **********************************************************************
6797 /* Reduction of degree when the squares are the squares of constraints. */
6801 /* TOUS,AB_SPECIFI::CARREAU&,REDUCTION,&CARREAU */
6803 /* INPUT ARGUMENTS : */
6804 /* ------------------ */
6805 /* NCFMXU: Max Nb of coeff by u of solution P(u,v) (table */
6806 /* PATCAN). This argument serves only to declare the size of this table. */
6807 /* NCFMXV: Max Nb of coeff by v of solution P(u,v) (table */
6808 /* PATCAN). This argument serves only to declare the size of this table. */
6809 /* NDIMEN: Total dimension of the space where the processed function */
6810 /* takes its values.(sum of dimensions of sub-spaces) */
6811 /* NBSESP: Nb of independent sub-spaces where the errors are measured. */
6812 /* NDIMSE: Table of dimensions of NBSESP sub-spaces. */
6813 /* NBUPAT: Nb of square solution by u. */
6814 /* NBVPAT: Nb of square solution by v. */
6815 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
6816 /* = 0, the extremities of iso-V are calculated */
6817 /* = 1, additionally the 1st derivative in the direction of iso-V is calculated */
6818 /* = 2, additionally the 2nd derivative in the direction of iso-V is calculated */
6819 /* IORDRV: Ordre de contrainte impose aux extremites de l'iso-U */
6820 /* = 0, on calcule les extremites de l'iso-U. */
6821 /* = 1, additionally the 1st derivative in the direction of iso-U is calculated */
6822 /* = 2, additionally the 2nd derivative in the direction of iso-U is calculated */
6823 /* EPSAPR: Table of imposed precisions, sub-space by sub-space. */
6824 /* EPSFRO: Table of imposed precisions, sub-space by sub-space on the limits of squares. */
6825 /* PATCAN: Table of coeff. in the canonic base of squares P(u,v) calculated for (u,v) in (-1,1). */
6826 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6827 /* committed in the approximation of F(u,v) by P(u,v). */
6828 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6829 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6831 /* OUTPUT ARGUMENTS : */
6832 /* ------------------- */
6833 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6834 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6836 /* COMMONS USED : */
6837 /* ---------------- */
6839 /* REFERENCES CALLED : */
6840 /* --------------------- */
6842 /* DESCRIPTION/NOTES/LIMITATIONS : */
6843 /* ------------------------------- */
6845 /* **********************************************************************
6848 /* Name of the routine */
6851 /* Parameter adjustments */
6852 epsfro_dim1 = *nbsesp;
6853 epsfro_offset = epsfro_dim1 * 5 + 1;
6854 epsfro -= epsfro_offset;
6857 ncoefv_dim1 = *nbupat;
6858 ncoefv_offset = ncoefv_dim1 + 1;
6859 ncoefv -= ncoefv_offset;
6860 ncoefu_dim1 = *nbupat;
6861 ncoefu_offset = ncoefu_dim1 + 1;
6862 ncoefu -= ncoefu_offset;
6863 errmax_dim1 = *nbsesp;
6864 errmax_dim2 = *nbupat;
6865 errmax_offset = errmax_dim1 * (errmax_dim2 + 1) + 1;
6866 errmax -= errmax_offset;
6867 patcan_dim1 = *ncfmxu;
6868 patcan_dim2 = *ncfmxv;
6869 patcan_dim3 = *ndimen;
6870 patcan_dim4 = *nbupat;
6871 patcan_offset = patcan_dim1 * (patcan_dim2 * (patcan_dim3 * (patcan_dim4
6873 patcan -= patcan_offset;
6876 ibb = AdvApp2Var_SysBase::mnfndeb_();
6878 AdvApp2Var_SysBase::mgenmsg_("MMA2FX6", 7L);
6883 for (jj = 1; jj <= i__1; ++jj) {
6885 for (ii = 1; ii <= i__2; ++ii) {
6886 ncfu = ncoefu[ii + jj * ncoefu_dim1];
6887 ncfv = ncoefv[ii + jj * ncoefv_dim1];
6889 /* ********************************************************************** */
6890 /* -------------------- Reduction of degree by U ------------------------- */
6891 /* ********************************************************************** */
6894 if (ncfu <= (*iordru + 1) << 1 && ncfu > 2) {
6898 for (ns = 1; ns <= i__3; ++ns) {
6901 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6902 tol = min(d__1,d__2);
6904 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6905 tol = min(d__1,d__2);
6907 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6908 tol = min(d__1,d__2);
6910 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6911 tol = min(d__1,d__2);
6912 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6915 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6916 tol = min(d__1,d__2);
6918 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6919 tol = min(d__1,d__2);
6921 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6922 tol = min(d__1,d__2);
6924 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6925 tol = min(d__1,d__2);
6930 for (nd = 1; nd <= i__4; ++nd) {
6933 for (kv = 1; kv <= i__5; ++kv) {
6934 bid += (d__1 = patcan[ncfu + (kv + (id + (ii + jj
6935 * patcan_dim4) * patcan_dim3) *
6936 patcan_dim2) * patcan_dim1], abs(d__1));
6942 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6943 errmax_dim2) * errmax_dim1]) {
6954 /* ********************************************************************** */
6955 /* -------------------- Reduction of degree by V ------------------------- */
6956 /* ********************************************************************** */
6959 if (ncfv <= (*iordrv + 1) << 1 && ncfv > 2) {
6963 for (ns = 1; ns <= i__3; ++ns) {
6966 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6967 tol = min(d__1,d__2);
6969 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6970 tol = min(d__1,d__2);
6972 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6973 tol = min(d__1,d__2);
6975 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6976 tol = min(d__1,d__2);
6977 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6980 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6981 tol = min(d__1,d__2);
6983 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6984 tol = min(d__1,d__2);
6986 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6987 tol = min(d__1,d__2);
6989 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6990 tol = min(d__1,d__2);
6995 for (nd = 1; nd <= i__4; ++nd) {
6998 for (ku = 1; ku <= i__5; ++ku) {
6999 bid += (d__1 = patcan[ku + (ncfv + (id + (ii + jj
7000 * patcan_dim4) * patcan_dim3) *
7001 patcan_dim2) * patcan_dim1], abs(d__1));
7007 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
7008 errmax_dim2) * errmax_dim1]) {
7019 /* --- Return the nbs of coeff. and pass to the next square --- */
7022 ncoefu[ii + jj * ncoefu_dim1] = max(ncfu,2);
7023 ncoefv[ii + jj * ncoefv_dim1] = max(ncfv,2);
7029 /* ------------------------------ The End -------------------------------
7033 AdvApp2Var_SysBase::mgsomsg_("MMA2FX6", 7L);
7039 //=======================================================================
7040 //function : mma2jmx_
7042 //=======================================================================
7043 int AdvApp2Var_ApproxF2var::mma2jmx_(integer *ndgjac,
7047 /* Initialized data */
7049 static doublereal xmax2[57] = { .9682458365518542212948163499456,
7050 .986013297183269340427888048593603,
7051 1.07810420343739860362585159028115,
7052 1.17325804490920057010925920756025,
7053 1.26476561266905634732910520370741,
7054 1.35169950227289626684434056681946,
7055 1.43424378958284137759129885012494,
7056 1.51281316274895465689402798226634,
7057 1.5878364329591908800533936587012,
7058 1.65970112228228167018443636171226,
7059 1.72874345388622461848433443013543,
7060 1.7952515611463877544077632304216,
7061 1.85947199025328260370244491818047,
7062 1.92161634324190018916351663207101,
7063 1.98186713586472025397859895825157,
7064 2.04038269834980146276967984252188,
7065 2.09730119173852573441223706382076,
7066 2.15274387655763462685970799663412,
7067 2.20681777186342079455059961912859,
7068 2.25961782459354604684402726624239,
7069 2.31122868752403808176824020121524,
7070 2.36172618435386566570998793688131,
7071 2.41117852396114589446497298177554,
7072 2.45964731268663657873849811095449,
7073 2.50718840313973523778244737914028,
7074 2.55385260994795361951813645784034,
7075 2.59968631659221867834697883938297,
7076 2.64473199258285846332860663371298,
7077 2.68902863641518586789566216064557,
7078 2.73261215675199397407027673053895,
7079 2.77551570192374483822124304745691,
7080 2.8177699459714315371037628127545,
7081 2.85940333797200948896046563785957,
7082 2.90044232019793636101516293333324,
7083 2.94091151970640874812265419871976,
7084 2.98083391718088702956696303389061,
7085 3.02023099621926980436221568258656,
7086 3.05912287574998661724731962377847,
7087 3.09752842783622025614245706196447,
7088 3.13546538278134559341444834866301,
7089 3.17295042316122606504398054547289,
7090 3.2099992681699613513775259670214,
7091 3.24662674946606137764916854570219,
7092 3.28284687953866689817670991319787,
7093 3.31867291347259485044591136879087,
7094 3.35411740487202127264475726990106,
7095 3.38919225660177218727305224515862,
7096 3.42390876691942143189170489271753,
7097 3.45827767149820230182596660024454,
7098 3.49230918177808483937957161007792,
7099 3.5260130200285724149540352829756,
7100 3.55939845146044235497103883695448,
7101 3.59247431368364585025958062194665,
7102 3.62524904377393592090180712976368,
7103 3.65773070318071087226169680450936,
7104 3.68992700068237648299565823810245,
7105 3.72184531357268220291630708234186 };
7106 static doublereal xmax4[55] = { 1.1092649593311780079813740546678,
7107 1.05299572648705464724876659688996,
7108 1.0949715351434178709281698645813,
7109 1.15078388379719068145021100764647,
7110 1.2094863084718701596278219811869,
7111 1.26806623151369531323304177532868,
7112 1.32549784426476978866302826176202,
7113 1.38142537365039019558329304432581,
7114 1.43575531950773585146867625840552,
7115 1.48850442653629641402403231015299,
7116 1.53973611681876234549146350844736,
7117 1.58953193485272191557448229046492,
7118 1.63797820416306624705258190017418,
7119 1.68515974143594899185621942934906,
7120 1.73115699602477936547107755854868,
7121 1.77604489805513552087086912113251,
7122 1.81989256661534438347398400420601,
7123 1.86276344480103110090865609776681,
7124 1.90471563564740808542244678597105,
7125 1.94580231994751044968731427898046,
7126 1.98607219357764450634552790950067,
7127 2.02556989246317857340333585562678,
7128 2.06433638992049685189059517340452,
7129 2.10240936014742726236706004607473,
7130 2.13982350649113222745523925190532,
7131 2.17661085564771614285379929798896,
7132 2.21280102016879766322589373557048,
7133 2.2484214321456956597803794333791,
7134 2.28349755104077956674135810027654,
7135 2.31805304852593774867640120860446,
7136 2.35210997297725685169643559615022,
7137 2.38568889602346315560143377261814,
7138 2.41880904328694215730192284109322,
7139 2.45148841120796359750021227795539,
7140 2.48374387161372199992570528025315,
7141 2.5155912654873773953959098501893,
7142 2.54704548720896557684101746505398,
7143 2.57812056037881628390134077704127,
7144 2.60882970619319538196517982945269,
7145 2.63918540521920497868347679257107,
7146 2.66919945330942891495458446613851,
7147 2.69888301230439621709803756505788,
7148 2.72824665609081486737132853370048,
7149 2.75730041251405791603760003778285,
7150 2.78605380158311346185098508516203,
7151 2.81451587035387403267676338931454,
7152 2.84269522483114290814009184272637,
7153 2.87060005919012917988363332454033,
7154 2.89823818258367657739520912946934,
7155 2.92561704377132528239806135133273,
7156 2.95274375377994262301217318010209,
7157 2.97962510678256471794289060402033,
7158 3.00626759936182712291041810228171,
7159 3.03267744830655121818899164295959,
7160 3.05886060707437081434964933864149 };
7161 static doublereal xmax6[53] = { 1.21091229812484768570102219548814,
7162 1.11626917091567929907256116528817,
7163 1.1327140810290884106278510474203,
7164 1.1679452722668028753522098022171,
7165 1.20910611986279066645602153641334,
7166 1.25228283758701572089625983127043,
7167 1.29591971597287895911380446311508,
7168 1.3393138157481884258308028584917,
7169 1.3821288728999671920677617491385,
7170 1.42420414683357356104823573391816,
7171 1.46546895108549501306970087318319,
7172 1.50590085198398789708599726315869,
7173 1.54550385142820987194251585145013,
7174 1.58429644271680300005206185490937,
7175 1.62230484071440103826322971668038,
7176 1.65955905239130512405565733793667,
7177 1.69609056468292429853775667485212,
7178 1.73193098017228915881592458573809,
7179 1.7671112206990325429863426635397,
7180 1.80166107681586964987277458875667,
7181 1.83560897003644959204940535551721,
7182 1.86898184653271388435058371983316,
7183 1.90180515174518670797686768515502,
7184 1.93410285411785808749237200054739,
7185 1.96589749778987993293150856865539,
7186 1.99721027139062501070081653790635,
7187 2.02806108474738744005306947877164,
7188 2.05846864831762572089033752595401,
7189 2.08845055210580131460156962214748,
7190 2.11802334209486194329576724042253,
7191 2.14720259305166593214642386780469,
7192 2.17600297710595096918495785742803,
7193 2.20443832785205516555772788192013,
7194 2.2325216999457379530416998244706,
7195 2.2602654243075083168599953074345,
7196 2.28768115912702794202525264301585,
7197 2.3147799369092684021274946755348,
7198 2.34157220782483457076721300512406,
7199 2.36806787963276257263034969490066,
7200 2.39427635443992520016789041085844,
7201 2.42020656255081863955040620243062,
7202 2.44586699364757383088888037359254,
7203 2.47126572552427660024678584642791,
7204 2.49641045058324178349347438430311,
7205 2.52130850028451113942299097584818,
7206 2.54596686772399937214920135190177,
7207 2.5703922285006754089328998222275,
7208 2.59459096001908861492582631591134,
7209 2.61856915936049852435394597597773,
7210 2.64233265984385295286445444361827,
7211 2.66588704638685848486056711408168,
7212 2.68923766976735295746679957665724,
7213 2.71238965987606292679677228666411 };
7215 /* System generated locals */
7218 /* Local variables */
7219 static logical ldbg;
7220 static integer numax, ii;
7221 static doublereal bid;
7224 /* **********************************************************************
7229 /* Calculate the max of Jacobo polynoms multiplied by the weight on */
7230 /* (-1,1) for order 0,4,6 or Legendre. */
7234 /* LEGENDRE,APPROXIMATION,ERREUR. */
7236 /* INPUT ARGUMENTS : */
7237 /* ------------------ */
7238 /* NDGJAC: Nb of Jacobi coeff. of approximation. */
7239 /* IORDRE: Order of continuity (from -1 to 2) */
7241 /* OUTPUT ARGUMENTS : */
7242 /* ------------------- */
7243 /* XJACMX: Table of maximums of Jacobi polynoms. */
7245 /* COMMONS USED : */
7246 /* ---------------- */
7248 /* REFERENCES CALLED : */
7249 /* --------------------- */
7251 /* DESCRIPTION/NOTES/LIMITATIONS : */
7252 /* ----------------------------------- */
7255 /* ***********************************************************************
7257 /* Name of the routine */
7258 /* ----------------------------- Initialisations ------------------------
7261 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7263 AdvApp2Var_SysBase::mgenmsg_("MMA2JMX", 7L);
7266 numax = *ndgjac - ((*iordre + 1) << 1);
7267 if (*iordre == -1) {
7269 for (ii = 0; ii <= i__1; ++ii) {
7270 bid = (ii * 2. + 1.) / 2.;
7271 xjacmx[ii] = sqrt(bid);
7274 } else if (*iordre == 0) {
7276 for (ii = 0; ii <= i__1; ++ii) {
7277 xjacmx[ii] = xmax2[ii];
7280 } else if (*iordre == 1) {
7282 for (ii = 0; ii <= i__1; ++ii) {
7283 xjacmx[ii] = xmax4[ii];
7286 } else if (*iordre == 2) {
7288 for (ii = 0; ii <= i__1; ++ii) {
7289 xjacmx[ii] = xmax6[ii];
7294 /* ------------------------- The end ------------------------------------
7298 AdvApp2Var_SysBase::mgsomsg_("MMA2JMX", 7L);
7303 //=======================================================================
7304 //function : mma2moy_
7306 //=======================================================================
7307 int mma2moy_(integer *ndgumx,
7319 /* System generated locals */
7320 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
7322 /* Local variables */
7323 static logical ldbg;
7324 static integer minu, minv, idebu, idebv, ii, nd, jj;
7325 static doublereal bid0, bid1;
7328 /* **********************************************************************
7333 /* Calculate the average approximation error made when only */
7334 /* the coefficients of PATJAC of degree between */
7335 /* 2*(IORDRU+1) and MINDGU by U and 2*(IORDRV+1) and MINDGV by V are preserved. */
7339 /* LEGENDRE,APPROXIMATION, AVERAGE ERROR */
7341 /* INPUT ARGUMENTS : */
7342 /* ------------------ */
7343 /* NDGUMX: Dimension by U of table PATJAC. */
7344 /* NDGVMX: Dimension by V of table PATJAC. */
7345 /* NDIMEN: Dimension of the space. */
7346 /* MINDGU: Lower limit of the index by U of PATJAC coeff to be taken into account. */
7347 /* MAXDGU: Upper limit of the index by U of PATJAC coeff to be taken into account. */
7348 /* MINDGV: Lower limit of the index by V of PATJAC coeff to be taken into account. */
7349 /* MAXDGV: Upper limit of the index by V of PATJAC coeff to be taken into account. */
7350 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
7351 /* IORDRV: Order of continuity by V provided by square PATJAC (from -1 to 2) */
7352 /* PATJAC: Table of coeff. of the approximation square with */
7353 /* constraints of order IORDRU by U and IORDRV by V. */
7355 /* OUTPUT ARGUMENTS : */
7356 /* ------------------- */
7357 /* ERRMOY: Average error commited by preserving only the coeff of */
7358 /* PATJAC 2*(IORDRU+1) in MINDGU by U and 2*(IORDRV+1) in MINDGV by V. */
7360 /* COMMONS USED : */
7361 /* ---------------- */
7363 /* REFERENCES CALLED : */
7364 /* --------------------- */
7366 /* DESCRIPTION/NOTES/LIMITATIONS : */
7367 /* ----------------------------------- */
7368 /* Table PATJAC stores the coeff. Cij of */
7369 /* approximation square F(U,V). Indexes i and j show the degree by */
7370 /* U and by V of the base polynoms. These base polynoms are in the form: */
7372 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
7374 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
7375 /* IORDRU+1 (the same by V by replacing U by V in the above expression). */
7377 /* The contribution to the average error of term Cij when */
7378 /* it is removed from PATJAC is Cij*Cij. */
7381 /* ***********************************************************************
7383 /* Name of the routine */
7386 /* ----------------------------- Initialisations ------------------------
7389 /* Parameter adjustments */
7390 patjac_dim1 = *ndgumx + 1;
7391 patjac_dim2 = *ndgvmx + 1;
7392 patjac_offset = patjac_dim1 * patjac_dim2;
7393 patjac -= patjac_offset;
7396 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7398 AdvApp2Var_SysBase::mgenmsg_("MMA2MOY", 7L);
7401 idebu = (*iordru + 1) << 1;
7402 idebv = (*iordrv + 1) << 1;
7403 minu = max(idebu,*mindgu);
7404 minv = max(idebv,*mindgv);
7408 /* ------------------ Calculation of the upper bound of the average error ------------ */
7409 /* -------------------- when the coeff. of indexes from MINDGU to MAXDGU ------ */
7410 /* ---------------- by U and of indexes from MINDGV to MAXDGV by V are removed -------------- */
7413 for (nd = 1; nd <= i__1; ++nd) {
7415 for (jj = minv; jj <= i__2; ++jj) {
7417 for (ii = idebu; ii <= i__3; ++ii) {
7418 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7419 bid0 += bid1 * bid1;
7428 for (nd = 1; nd <= i__1; ++nd) {
7430 for (jj = idebv; jj <= i__2; ++jj) {
7432 for (ii = minu; ii <= i__3; ++ii) {
7433 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7434 bid0 += bid1 * bid1;
7442 /* ----------------------- Calculation of the average error -------------
7446 *errmoy = sqrt(bid0);
7448 /* ------------------------- The end ------------------------------------
7452 AdvApp2Var_SysBase::mgsomsg_("MMA2MOY", 7L);
7457 //=======================================================================
7458 //function : mma2roo_
7460 //=======================================================================
7461 int AdvApp2Var_ApproxF2var::mma2roo_(integer *nbpntu,
7466 /* System generated locals */
7469 /* Local variables */
7470 static integer ii, ibb;
7472 /* **********************************************************************
7477 /* Return roots of Legendre for discretisations. */
7481 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
7483 /* INPUT ARGUMENTS : */
7484 /* ------------------ */
7485 /* NBPNTU: Nb of INTERNAL parameters of discretization BY U. */
7486 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7487 /* NBPNTV: Nb of INTERNAL parameters of discretization BY V. */
7488 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7490 /* OUTPUT ARGUMENTS : */
7491 /* ------------------- */
7492 /* UROOTL: Table of parameters of discretisation ON (-1,1) BY U.
7494 /* VROOTL: Table of parameters of discretisation ON (-1,1) BY V.
7497 /* COMMONS USED : */
7498 /* ---------------- */
7500 /* REFERENCES CALLED : */
7501 /* --------------------- */
7503 /* DESCRIPTION/NOTES/LIMITATIONS : */
7504 /* ----------------------------------- */
7507 /* **********************************************************************
7510 /* Name of the routine */
7513 /* Parameter adjustments */
7518 ibb = AdvApp2Var_SysBase::mnfndeb_();
7520 AdvApp2Var_SysBase::mgenmsg_("MMA2ROO", 7L);
7523 /* ---------------- Return the POSITIVE roots on U ------------------
7526 AdvApp2Var_MathBase::mmrtptt_(nbpntu, &urootl[(*nbpntu + 1) / 2 + 1]);
7528 for (ii = 1; ii <= i__1; ++ii) {
7529 urootl[ii] = -urootl[*nbpntu - ii + 1];
7532 if (*nbpntu % 2 == 1) {
7533 urootl[*nbpntu / 2 + 1] = 0.;
7536 /* ---------------- Return the POSITIVE roots on V ------------------
7539 AdvApp2Var_MathBase::mmrtptt_(nbpntv, &vrootl[(*nbpntv + 1) / 2 + 1]);
7541 for (ii = 1; ii <= i__1; ++ii) {
7542 vrootl[ii] = -vrootl[*nbpntv - ii + 1];
7545 if (*nbpntv % 2 == 1) {
7546 vrootl[*nbpntv / 2 + 1] = 0.;
7549 /* ------------------------------ The End -------------------------------
7553 AdvApp2Var_SysBase::mgsomsg_("MMA2ROO", 7L);
7557 //=======================================================================
7558 //function : mmmapcoe_
7560 //=======================================================================
7561 int mmmapcoe_(integer *ndim,
7571 /* System generated locals */
7572 integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
7573 crvjac_dim1, crvjac_offset, gsstab_dim1, i__1, i__2, i__3;
7575 /* Local variables */
7576 static integer igss, ikdeb;
7577 static doublereal bidon;
7578 static integer nd, ik, ir, nbroot, ibb;
7582 /* **********************************************************************
7587 /* Calculate the coefficients of polinomial approximation curve */
7588 /* of degree NDGJAC by the method of smallest squares starting from */
7589 /* the discretization of function on the roots of Legendre polynom */
7590 /* of degree NBPNTS. */
7594 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
7596 /* INPUT ARGUMENTS : */
7597 /* ------------------ */
7598 /* NDIM : Dimension of the space. */
7599 /* NDGJAC : Max Degree of the polynom of approximation. */
7600 /* The representation in the orthogonal base starts from degree */
7601 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7602 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7603 /* IORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7604 /* to step of constraints, C0,C1 or C2. */
7605 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7606 /* are calculated the coefficients of integration by */
7607 /* Gauss method. It is required to set NBPNTS=30,40,50 or 61 */
7608 /* and NDGJAC < NBPNTS. */
7609 /* SOMTAB : Table of F(ti)+F(-ti) with ti in ROOTAB. */
7610 /* DIFTAB : Table of F(ti)-F(-ti) with ti in ROOTAB. */
7611 /* GSSTAB(i,k) : Table of coefficients of integration by the Gauss method : */
7612 /* i varies from 0 to NBPNTS and */
7613 /* k varies from 0 to NDGJAC-2*(JORDRE+1). */
7615 /* OUTPUT ARGUMENTSE : */
7616 /* ------------------- */
7617 /* CRVJAC : Curve of approximation of FONCNP with eventually */
7618 /* taking into account of constraints at the extremities. */
7619 /* This curve is of degree NDGJAC. */
7621 /* COMMONS USED : */
7622 /* ---------------- */
7624 /* REFERENCES CALLED : */
7625 /* --------------------- */
7627 /* DESCRIPTION/NOTES/LIMITATIONS : */
7628 /* ------------------------------- */
7630 /* **********************************************************************
7633 /* Name of the routine */
7635 /* Parameter adjustments */
7636 crvjac_dim1 = *ndgjac + 1;
7637 crvjac_offset = crvjac_dim1;
7638 crvjac -= crvjac_offset;
7639 gsstab_dim1 = *nbpnts / 2 + 1;
7640 diftab_dim1 = *nbpnts / 2 + 1;
7641 diftab_offset = diftab_dim1;
7642 diftab -= diftab_offset;
7643 somtab_dim1 = *nbpnts / 2 + 1;
7644 somtab_offset = somtab_dim1;
7645 somtab -= somtab_offset;
7648 ibb = AdvApp2Var_SysBase::mnfndeb_();
7650 AdvApp2Var_SysBase::mgenmsg_("MMMAPCO", 7L);
7652 ikdeb = (*iordre + 1) << 1;
7653 nbroot = *nbpnts / 2;
7656 for (nd = 1; nd <= i__1; ++nd) {
7658 /* ----------------- Calculate the coefficients of even degree ----------
7662 for (ik = ikdeb; ik <= i__2; ik += 2) {
7666 for (ir = 1; ir <= i__3; ++ir) {
7667 bidon += somtab[ir + nd * somtab_dim1] * gsstab[ir + igss *
7671 crvjac[ik + nd * crvjac_dim1] = bidon;
7675 /* --------------- Calculate the coefficients of uneven degree ----------
7679 for (ik = ikdeb + 1; ik <= i__2; ik += 2) {
7683 for (ir = 1; ir <= i__3; ++ir) {
7684 bidon += diftab[ir + nd * diftab_dim1] * gsstab[ir + igss *
7688 crvjac[ik + nd * crvjac_dim1] = bidon;
7695 /* ------- Add terms connected to the supplementary root (0.D0) ------ */
7696 /* ----------- of Legendre polynom of uneven degree NBPNTS -----------
7699 if (*nbpnts % 2 == 0) {
7703 for (nd = 1; nd <= i__1; ++nd) {
7705 for (ik = ikdeb; ik <= i__2; ik += 2) {
7707 crvjac[ik + nd * crvjac_dim1] += somtab[nd * somtab_dim1] *
7708 gsstab[igss * gsstab_dim1];
7714 /* ------------------------------ The end -------------------------------
7719 AdvApp2Var_SysBase::mgsomsg_("MMMAPCO", 7L);
7723 //=======================================================================
7724 //function : mmaperm_
7726 //=======================================================================
7727 int mmaperm_(integer *ncofmx,
7735 /* System generated locals */
7736 integer crvjac_dim1, crvjac_offset, i__1, i__2;
7738 /* Local variables */
7739 static doublereal bidj;
7740 static integer i__, ia, nd, ncfcut, ibb;
7741 static doublereal bid;
7745 /* **********************************************************************
7750 /* Calculate the square root of the average quadratic error */
7751 /* of approximation done when only the */
7752 /* first NCFNEW coefficients of a curve of degree NCOEFF-1 */
7753 /* written in NORMALIZED Jacobi base of order 2*(IORDRE+1) are preserved. */
7757 /* LEGENDRE,POLYGONE,APPROXIMATION,ERREUR. */
7759 /* INPUT ARGUMENTS : */
7760 /* ------------------ */
7761 /* NCOFMX : Maximum degree of the curve. */
7762 /* NDIM : Dimension of the space. */
7763 /* NCOEFF : Degree +1 of the curve. */
7764 /* IORDRE : Order of constraint of continuity at the extremities. */
7765 /* CRVJAC : The curve the degree which of will be lowered. */
7766 /* NCFNEW : Degree +1 of the resulting polynom. */
7768 /* OUTPUT ARGUMENTS : */
7769 /* ------------------- */
7770 /* ERRMOY : Average precision of approximation. */
7772 /* COMMONS USED : */
7773 /* ---------------- */
7775 /* REFERENCES CALLED : */
7776 /* ----------------------- */
7778 /* DESCRIPTION/NOTES/LIMITATIONS : */
7779 /* ----------------------------------- */
7781 /* ***********************************************************************
7784 /* Name of the routine */
7786 /* Parameter adjustments */
7787 crvjac_dim1 = *ncofmx;
7788 crvjac_offset = crvjac_dim1 + 1;
7789 crvjac -= crvjac_offset;
7792 ibb = AdvApp2Var_SysBase::mnfndeb_();
7794 AdvApp2Var_SysBase::mgenmsg_("MMAPERM", 7L);
7797 /* --------- Minimum degree that can be reached : Stop at 1 or IA -------
7800 ia = (*iordre + 1) << 1;
7802 if (*ncfnew + 1 > ncfcut) {
7803 ncfcut = *ncfnew + 1;
7806 /* -------------- Elimination of coefficients of high degree ------------ */
7807 /* ----------- Loop on the series of Jacobi :NCFCUT --> NCOEFF --------- */
7812 for (nd = 1; nd <= i__1; ++nd) {
7814 for (i__ = ncfcut; i__ <= i__2; ++i__) {
7815 bidj = crvjac[i__ + nd * crvjac_dim1];
7822 /* ----------- Square Root of average quadratic error e -----------
7826 *errmoy = sqrt(bid);
7828 /* ------------------------------- The end ------------------------------
7832 AdvApp2Var_SysBase::mgsomsg_("MMAPERM", 7L);
7836 //=======================================================================
7837 //function : mmapptt_
7839 //=======================================================================
7840 int AdvApp2Var_ApproxF2var::mmapptt_(const integer *ndgjac,
7841 const integer *nbpnts,
7842 const integer *jordre,
7846 /* System generated locals */
7847 integer cgauss_dim1, i__1;
7849 /* Local variables */
7850 static integer kjac, iptt, ipdb0, infdg, iptdb, mxjac, ilong, ibb;
7854 /* **********************************************************************
7859 /* Load the elements required for integration by */
7860 /* Gauss method to obtain the coefficients in the base of
7861 /* Legendre of the approximation by the least squares of a */
7862 /* function. The elements are stored in commons MMAPGSS */
7863 /* (case without constraint), MMAPGS0 (constraints C0), MMAPGS1 */
7864 /* (constraints C1) and MMAPGS2 (constraints C2). */
7868 /* INTEGRATION,GAUSS,JACOBI */
7870 /* INPUT ARGUMENTS : */
7871 /* ------------------ */
7872 /* NDGJAC : Max degree of the polynom of approximation. */
7873 /* The representation in orthogonal base goes from degree
7874 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7875 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7876 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7877 /* are calculated the coefficients of integration by the */
7878 /* method of Gauss. It is required that NBPNTS=8,10,15,20,25, */
7879 /* 30,40,50 or 61 and NDGJAC < NBPNTS. */
7880 /* JORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7881 /* to step of constraints C0,C1 or C2. */
7883 /* OUTPUT ARGUMENTS : */
7884 /* ------------------- */
7885 /* CGAUSS(i,k) : Table of coefficients of integration by */
7886 /* Gauss method : i varies from 0 to the integer part */
7887 /* of NBPNTS/2 and k varies from 0 to NDGJAC-2*(JORDRE+1). */
7888 /* These are the coeff. of integration associated to */
7889 /* positive roots of the polynom of Legendre of degree */
7890 /* NBPNTS. CGAUSS(0,k) contains coeff. */
7891 /* of integration associated to root t = 0 when */
7892 /* NBPNTS is uneven. */
7893 /* IERCOD : Error code. */
7895 /* = 11 NBPNTS is not 8,10,15,20,25,30,40,50 or 61. */
7896 /* = 21 JORDRE is not -1,0,1 or 2. */
7897 /* = 31 NDGJAC is too great or too small. */
7899 /* COMMONS USED : */
7900 /* ---------------- */
7901 /* MMAPGSS,MMAPGS0,MMAPGS1,MMAPGS2. */
7902 /* ***********************************************************************
7904 /* Parameter adjustments */
7905 cgauss_dim1 = *nbpnts / 2 + 1;
7908 ibb = AdvApp2Var_SysBase::mnfndeb_();
7910 AdvApp2Var_SysBase::mgenmsg_("MMAPPTT", 7L);
7914 /* ------------------- Tests on the validity of inputs ----------------
7917 infdg = (*jordre + 1) << 1;
7918 if (*nbpnts != 8 && *nbpnts != 10 && *nbpnts != 15 && *nbpnts != 20 && *
7919 nbpnts != 25 && *nbpnts != 30 && *nbpnts != 40 && *nbpnts != 50 &&
7924 if (*jordre < -1 || *jordre > 2) {
7928 if (*ndgjac >= *nbpnts || *ndgjac < infdg) {
7932 /* --------------- Calculation of the start pointer following NBPNTS -----------
7937 iptdb += (8 - infdg) << 2;
7940 iptdb += (10 - infdg) * 5;
7943 iptdb += (15 - infdg) * 7;
7946 iptdb += (20 - infdg) * 10;
7949 iptdb += (25 - infdg) * 12;
7952 iptdb += (30 - infdg) * 15;
7955 iptdb += (40 - infdg) * 20;
7958 iptdb += (50 - infdg) * 25;
7963 ipdb0 = ipdb0 + (14 - infdg) / 2 + 1;
7966 ipdb0 = ipdb0 + (24 - infdg) / 2 + 1;
7969 /* ------------------ Choice of the common depending on JORDRE -------------
7972 if (*jordre == -1) {
7985 /* ---------------- Common MMAPGSS (case without constraints) ----------------
7989 ilong = *nbpnts / 2 << 3;
7991 for (kjac = 0; kjac <= i__1; ++kjac) {
7992 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7993 AdvApp2Var_SysBase::mcrfill_(&ilong,
7994 (char *)&mmapgss_.gslxjs[iptt - 1],
7995 (char *)&cgauss[kjac * cgauss_dim1 + 1]);
7998 /* --> Case when the number of points is uneven. */
7999 if (*nbpnts % 2 == 1) {
8002 for (kjac = 0; kjac <= i__1; kjac += 2) {
8003 cgauss[kjac * cgauss_dim1] = mmapgss_.gsl0js[iptt - 1];
8008 for (kjac = 1; kjac <= i__1; kjac += 2) {
8009 cgauss[kjac * cgauss_dim1] = 0.;
8015 /* ---------------- Common MMAPGS0 (case with constraints C0) -------------
8019 mxjac = *ndgjac - infdg;
8020 ilong = *nbpnts / 2 << 3;
8022 for (kjac = 0; kjac <= i__1; ++kjac) {
8023 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
8024 AdvApp2Var_SysBase::mcrfill_(&ilong,
8025 (char *)&mmapgs0_.gslxj0[iptt - 1],
8026 (char *)&cgauss[kjac * cgauss_dim1 + 1]);
8029 /* --> Case when the number of points is uneven. */
8030 if (*nbpnts % 2 == 1) {
8033 for (kjac = 0; kjac <= i__1; kjac += 2) {
8034 cgauss[kjac * cgauss_dim1] = mmapgs0_.gsl0j0[iptt - 1];
8039 for (kjac = 1; kjac <= i__1; kjac += 2) {
8040 cgauss[kjac * cgauss_dim1] = 0.;
8046 /* ---------------- Common MMAPGS1 (case with constraints C1) -------------
8050 mxjac = *ndgjac - infdg;
8051 ilong = *nbpnts / 2 << 3;
8053 for (kjac = 0; kjac <= i__1; ++kjac) {
8054 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
8055 AdvApp2Var_SysBase::mcrfill_(&ilong,
8056 (char *)&mmapgs1_.gslxj1[iptt - 1],
8057 (char *)&cgauss[kjac * cgauss_dim1 + 1]);
8060 /* --> Case when the number of points is uneven. */
8061 if (*nbpnts % 2 == 1) {
8064 for (kjac = 0; kjac <= i__1; kjac += 2) {
8065 cgauss[kjac * cgauss_dim1] = mmapgs1_.gsl0j1[iptt - 1];
8070 for (kjac = 1; kjac <= i__1; kjac += 2) {
8071 cgauss[kjac * cgauss_dim1] = 0.;
8077 /* ---------------- Common MMAPGS2 (case with constraints C2) -------------
8081 mxjac = *ndgjac - infdg;
8082 ilong = *nbpnts / 2 << 3;
8084 for (kjac = 0; kjac <= i__1; ++kjac) {
8085 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
8086 AdvApp2Var_SysBase::mcrfill_(&ilong,
8087 (char *)&mmapgs2_.gslxj2[iptt - 1],
8088 (char *)&cgauss[kjac * cgauss_dim1 + 1]);
8091 /* --> Cas of uneven number of points. */
8092 if (*nbpnts % 2 == 1) {
8095 for (kjac = 0; kjac <= i__1; kjac += 2) {
8096 cgauss[kjac * cgauss_dim1] = mmapgs2_.gsl0j2[iptt - 1];
8101 for (kjac = 1; kjac <= i__1; kjac += 2) {
8102 cgauss[kjac * cgauss_dim1] = 0.;
8108 /* ------------------------- Return the error code --------------
8110 /* --> NBPNTS is not OK */
8114 /* --> JORDRE is not OK */
8118 /* --> NDGJAC is not OK */
8123 /* -------------------------------- The end -----------------------------
8128 AdvApp2Var_SysBase::maermsg_("MMAPPTT", iercod, 7L);
8131 AdvApp2Var_SysBase::mgsomsg_("MMAPPTT", 7L);
8137 //=======================================================================
8138 //function : mmjacpt_
8140 //=======================================================================
8141 int mmjacpt_(const integer *ndimen,
8142 const integer *ncoefu,
8143 const integer *ncoefv,
8144 const integer *iordru,
8145 const integer *iordrv,
8146 const doublereal *ptclgd,
8150 /* System generated locals */
8151 integer ptccan_dim1, ptccan_dim2, ptccan_offset, ptclgd_dim1, ptclgd_dim2,
8152 ptclgd_offset, ptcaux_dim1, ptcaux_dim2, ptcaux_dim3,
8153 ptcaux_offset, i__1, i__2, i__3;
8155 /* Local variables */
8156 static integer kdim, nd, ii, jj, ibb;
8158 /* ***********************************************************************
8163 /* Passage from canonical to Jacobi base for a */
8164 /* "square" in a space of arbitrary dimension. */
8168 /* SMOOTHING,BASE,LEGENDRE */
8171 /* INPUT ARGUMENTS : */
8172 /* ------------------ */
8173 /* NDIMEN : Dimension of the space. */
8174 /* NCOEFU : Degree+1 by U. */
8175 /* NCOEFV : Degree+1 by V. */
8176 /* IORDRU : Order of Jacobi polynoms by U. */
8177 /* IORDRV : Order of Jacobi polynoms by V. */
8178 /* PTCLGD : The square in the Jacobi base. */
8180 /* OUTPUT ARGUMENTS : */
8181 /* ------------------- */
8182 /* PTCAUX : Auxilliary space. */
8183 /* PTCCAN : The square in the canonic base (-1,1) */
8185 /* COMMONS USED : */
8186 /* ---------------- */
8188 /* APPLIED REFERENCES : */
8189 /* ----------------------- */
8191 /* DESCRIPTION/NOTES/LIMITATIONS : */
8192 /* ----------------------------------- */
8193 /* Cancels and replaces MJACPC */
8195 /* *********************************************************************
8197 /* Name of the routine */
8200 /* Parameter adjustments */
8201 ptccan_dim1 = *ncoefu;
8202 ptccan_dim2 = *ncoefv;
8203 ptccan_offset = ptccan_dim1 * (ptccan_dim2 + 1) + 1;
8204 ptccan -= ptccan_offset;
8205 ptcaux_dim1 = *ncoefv;
8206 ptcaux_dim2 = *ncoefu;
8207 ptcaux_dim3 = *ndimen;
8208 ptcaux_offset = ptcaux_dim1 * (ptcaux_dim2 * (ptcaux_dim3 + 1) + 1) + 1;
8209 ptcaux -= ptcaux_offset;
8210 ptclgd_dim1 = *ncoefu;
8211 ptclgd_dim2 = *ncoefv;
8212 ptclgd_offset = ptclgd_dim1 * (ptclgd_dim2 + 1) + 1;
8213 ptclgd -= ptclgd_offset;
8216 ibb = AdvApp2Var_SysBase::mnfndeb_();
8218 AdvApp2Var_SysBase::mgenmsg_("MMJACPT", 7L);
8221 /* Passage into canonical by u. */
8223 kdim = *ndimen * *ncoefv;
8224 AdvApp2Var_MathBase::mmjaccv_((integer *)ncoefu,
8227 (doublereal *)&ptclgd[ptclgd_offset],
8228 (doublereal *)&ptcaux[ptcaux_offset],
8229 (doublereal *)&ptccan[ptccan_offset]);
8231 /* Swapping of u and v. */
8234 for (nd = 1; nd <= i__1; ++nd) {
8236 for (jj = 1; jj <= i__2; ++jj) {
8238 for (ii = 1; ii <= i__3; ++ii) {
8239 ptcaux[jj + (ii + (nd + ptcaux_dim3) * ptcaux_dim2) *
8240 ptcaux_dim1] = ptccan[ii + (jj + nd * ptccan_dim2) *
8249 /* Passage into canonical by v. */
8251 kdim = *ndimen * *ncoefu;
8252 AdvApp2Var_MathBase::mmjaccv_((integer *)ncoefv,
8255 (doublereal *)&ptcaux[((ptcaux_dim3 + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1],
8256 (doublereal *)&ptccan[ptccan_offset],
8257 (doublereal *)&ptcaux[(((ptcaux_dim3 << 1) + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1]);
8259 /* Swapping of u and v. */
8262 for (nd = 1; nd <= i__1; ++nd) {
8264 for (jj = 1; jj <= i__2; ++jj) {
8266 for (ii = 1; ii <= i__3; ++ii) {
8267 ptccan[ii + (jj + nd * ptccan_dim2) * ptccan_dim1] = ptcaux[
8268 jj + (ii + (nd + (ptcaux_dim3 << 1)) * ptcaux_dim2) *
8277 /* ---------------------------- THAT'S ALL FOLKS ------------------------
8281 AdvApp2Var_SysBase::mgsomsg_("MMJACPT", 7L);