1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 // This file is part of Open CASCADE Technology software library.
5 // This library is free software; you can redistribute it and / or modify it
6 // under the terms of the GNU Lesser General Public version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
14 // AdvApp2Var_ApproxF2var.cxx
16 #include <AdvApp2Var_SysBase.hxx>
17 #include <AdvApp2Var_MathBase.hxx>
18 #include <AdvApp2Var_Data_f2c.hxx>
19 #include <AdvApp2Var_Data.hxx>
20 #include <AdvApp2Var_ApproxF2var.hxx>
24 int mmjacpt_(const integer *ndimen,
25 const integer *ncoefu,
26 const integer *ncoefv,
27 const integer *iordru,
28 const integer *iordrv,
29 const doublereal *ptclgd,
36 int mma2ce2_(integer *numdec,
71 int mma2cfu_(integer *ndujac,
83 int mma2cfv_(integer *ndvjac,
93 int mma2er1_(integer *ndjacu,
109 int mma2er2_(integer *ndjacu,
128 int mma2moy_(integer *ndgumx,
141 int mma2ds2_(integer *ndimen,
144 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
162 int mma1fdi_(integer *ndimen,
164 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
179 int mma1cdi_(integer *ndimen,
191 int mma1jak_(integer *ndimen,
201 int mma1cnt_(integer *ndimen,
210 int mma1fer_(integer *ndimen,
225 int mma1noc_(doublereal *dfuvin,
236 int mmmapcoe_(integer *ndim,
246 int mmaperm_(integer *ncofmx,
255 #define mmapgss_1 mmapgss_
256 #define mmapgs0_1 mmapgs0_
257 #define mmapgs1_1 mmapgs1_
258 #define mmapgs2_1 mmapgs2_
260 //=======================================================================
261 //function : mma1cdi_
263 //=======================================================================
264 int mma1cdi_(integer *ndimen,
278 /* System generated locals */
279 integer contr1_dim1, contr1_offset, contr2_dim1, contr2_offset,
280 somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
281 fpntab_dim1, fpntab_offset, hermit_dim1, hermit_offset, i__1,
284 /* Local variables */
285 integer nroo2, ncfhe, nd, ii, kk;
286 integer ibb, kkm, kkp;
287 doublereal bid1, bid2, bid3 = 0.;
289 /* **********************************************************************
293 /* Discretisation on the parameters of interpolation polynomes */
294 /* constraints of order IORDRE. */
298 /* ALL, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
300 /* INPUT ARGUMENTS : */
301 /* ------------------ */
302 /* NDIMEN: Space dimension. */
303 /* NBROOT: Number of INTERNAL discretisation parameters. */
304 /* It is also the root number Legendre polynome where */
305 /* the discretization is performed. */
306 /* ROOTLG: Table of discretization parameters ON (-1,1). */
307 /* IORDRE: Order of constraint imposed to the extremities of the iso. */
308 /* = 0, the extremities of the iso are calculated */
309 /* = 1, additionally, the 1st derivative in the direction */
310 /* of the iso is calculated. */
311 /* = 2, additionally, the 2nd derivative in the direction */
312 /* of the iso is calculated. */
313 /* CONTR1: Contains, if IORDRE>=0, values IORDRE+1 in TTABLE(0)
315 /* (1st extremity) of derivatives of F(Uc,Ve) or F(Ue,Vc), */
317 /* CONTR2: Contains, if IORDRE>=0, values IORDRE+1 in */
318 /* TTABLE(NBROOT+1) (2nd extremity) of: */
319 /* If ISOFAV=1, derived of order IDERIV by U, derived */
320 /* ordre 0 to IORDRE by V of F(Uc,Ve) or Uc=TCONST */
321 /* (fixed iso value) and Ve is the fixed extremity. */
322 /* If ISOFAV=2, derivative of order IDERIV by V, derivative */
323 /* of order 0 to IORDRE by U of F(Ue,Vc) or Vc=TCONST */
324 /* (fixed iso value) and Ue is the fixed extremity. */
326 /* SOMTAB: Table of NBROOT/2 sums of 2 index points */
327 /* NBROOT-II+1 and II, for II = 1, NBROOT/2. */
328 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
329 /* NBROOT-II+1 and II, for II = 1, NBROOT/2. */
331 /* OUTPUT ARGUMENTS : */
332 /* ------------------- */
333 /* SOMTAB: Table of NBROOT/2 sums of 2 index points */
334 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
335 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
336 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
337 /* FPNTAB: Auxiliary table. */
338 /* HERMIT: Table of coeff. 2*(IORDRE+1) Hermite polynoms */
339 /* of degree 2*IORDRE+1. */
340 /* IERCOD: Error code, */
341 /* = 0, Everythig is OK */
342 /* = 1, The value of IORDRE is out of (0,2) */
344 /* ---------------- */
346 /* REFERENCES CALLED : */
347 /* ----------------------- */
349 /* DESCRIPTION/NOTES/LIMITATIONS : */
350 /* ----------------------------------- */
351 /* The results of discretization are arranged in 2 tables */
352 /* SOMTAB and DIFTAB to earn time during the */
353 /* calculation of coefficients of the approximation curve. */
355 /* If NBROOT is uneven in SOMTAB(0,*) and DIFTAB(0,*) one stores */
356 /* the values of the median root of Legendre (0.D0 in (-1,1)). */
358 /* **********************************************************************
361 /* Name of the routine */
364 /* Parameter adjustments */
365 diftab_dim1 = *nbroot / 2 + 1;
366 diftab_offset = diftab_dim1;
367 diftab -= diftab_offset;
368 somtab_dim1 = *nbroot / 2 + 1;
369 somtab_offset = somtab_dim1;
370 somtab -= somtab_offset;
372 hermit_dim1 = (*iordre << 1) + 2;
373 hermit_offset = hermit_dim1;
374 hermit -= hermit_offset;
375 fpntab_dim1 = *nbroot;
376 fpntab_offset = fpntab_dim1 + 1;
377 fpntab -= fpntab_offset;
378 contr2_dim1 = *ndimen;
379 contr2_offset = contr2_dim1 + 1;
380 contr2 -= contr2_offset;
381 contr1_dim1 = *ndimen;
382 contr1_offset = contr1_dim1 + 1;
383 contr1 -= contr1_offset;
386 ibb = AdvApp2Var_SysBase::mnfndeb_();
388 AdvApp2Var_SysBase::mgenmsg_("MMA1CDI", 7L);
392 /* --- Recuperate 2*(IORDRE+1) coeff of 2*(IORDRE+1) of Hermite polynom ---
395 AdvApp2Var_ApproxF2var::mma1her_(iordre, &hermit[hermit_offset], iercod);
400 /* ------------------- Discretization of Hermite polynoms -----------
403 ncfhe = (*iordre + 1) << 1;
405 for (ii = 1; ii <= i__1; ++ii) {
407 for (kk = 1; kk <= i__2; ++kk) {
408 AdvApp2Var_MathBase::mmmpocur_(&ncfhe, &c__1, &ncfhe, &hermit[ii * hermit_dim1], &
409 rootlg[kk], &fpntab[kk + ii * fpntab_dim1]);
415 /* ---- Discretizations of boundary polynoms are taken ----
420 for (nd = 1; nd <= i__1; ++nd) {
422 for (ii = 1; ii <= i__2; ++ii) {
423 bid1 = contr1[nd + ii * contr1_dim1];
424 bid2 = contr2[nd + ii * contr2_dim1];
426 for (kk = 1; kk <= i__3; ++kk) {
427 kkm = nroo2 - kk + 1;
428 bid3 = bid1 * fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1] +
429 bid2 * fpntab[kkm + (ii << 1) * fpntab_dim1];
430 somtab[kk + nd * somtab_dim1] -= bid3;
431 diftab[kk + nd * diftab_dim1] += bid3;
435 for (kk = 1; kk <= i__3; ++kk) {
436 kkp = (*nbroot + 1) / 2 + kk;
437 bid3 = bid1 * fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] +
438 bid2 * fpntab[kkp + (ii << 1) * fpntab_dim1];
439 somtab[kk + nd * somtab_dim1] -= bid3;
440 diftab[kk + nd * diftab_dim1] -= bid3;
448 /* ------------ Cas when discretization is done on the roots of a -----------
450 /* ---------- Legendre polynom of uneven degree, 0 is root --------
453 if (*nbroot % 2 == 1) {
455 for (nd = 1; nd <= i__1; ++nd) {
457 for (ii = 1; ii <= i__2; ++ii) {
458 bid3 = fpntab[nroo2 + 1 + ((ii << 1) - 1) * fpntab_dim1] *
459 contr1[nd + ii * contr1_dim1] + fpntab[nroo2 + 1 + (
460 ii << 1) * fpntab_dim1] * contr2[nd + ii *
464 somtab[nd * somtab_dim1] -= bid3;
465 diftab[nd * diftab_dim1] -= bid3;
472 /* ------------------------------ The End -------------------------------
474 /* --> IORDRE is not in the authorized zone. */
481 AdvApp2Var_SysBase::mgsomsg_("MMA1CDI", 7L);
486 //=======================================================================
487 //function : mma1cnt_
489 //=======================================================================
490 int mma1cnt_(integer *ndimen,
498 /* System generated locals */
499 integer contr1_dim1, contr1_offset, contr2_dim1, contr2_offset,
500 hermit_dim1, hermit_offset, crvjac_dim1, crvjac_offset, i__1,
503 /* Local variables */
504 integer nd, ii, jj, ibb;
508 /* ***********************************************************************
513 /* Add constraint to polynom. */
517 /* ALL,AB_SPECIFI::COURE&,APPROXIMATION,ADDITION,&CONSTRAINT */
519 /* INPUT ARGUMENTS : */
520 /* -------------------- */
521 /* NDIMEN: Dimension of the space */
522 /* IORDRE: Order of constraint. */
523 /* CONTR1: pt of constraint in -1, from order 0 to IORDRE. */
524 /* CONTR2: Pt of constraint in +1, from order 0 to IORDRE. */
525 /* HERMIT: Table of Hermit polynoms of order IORDRE. */
526 /* CRVJAV: Curve of approximation in Jacobi base. */
528 /* OUTPUT ARGUMENTS : */
529 /* --------------------- */
530 /* CRVJAV: Curve of approximation in Jacobi base */
531 /* to which the polynom of interpolation of constraints is added. */
534 /* ------------------ */
537 /* REFERENCES CALLED : */
538 /* --------------------- */
541 /* DESCRIPTION/NOTES/LIMITATIONS : */
542 /* ----------------------------------- */
545 /* ***********************************************************************
548 /* ***********************************************************************
550 /* Name of the routine */
552 /* ***********************************************************************
554 /* INITIALISATIONS */
555 /* ***********************************************************************
558 /* Parameter adjustments */
559 hermit_dim1 = (*iordre << 1) + 2;
560 hermit_offset = hermit_dim1;
561 hermit -= hermit_offset;
562 contr2_dim1 = *ndimen;
563 contr2_offset = contr2_dim1 + 1;
564 contr2 -= contr2_offset;
565 contr1_dim1 = *ndimen;
566 contr1_offset = contr1_dim1 + 1;
567 contr1 -= contr1_offset;
568 crvjac_dim1 = *ndgjac + 1;
569 crvjac_offset = crvjac_dim1;
570 crvjac -= crvjac_offset;
573 ibb = AdvApp2Var_SysBase::mnfndeb_();
575 AdvApp2Var_SysBase::mgenmsg_("MMA1CNT", 7L);
578 /* ***********************************************************************
581 /* ***********************************************************************
585 for (nd = 1; nd <= i__1; ++nd) {
586 i__2 = (*iordre << 1) + 1;
587 for (ii = 0; ii <= i__2; ++ii) {
590 for (jj = 1; jj <= i__3; ++jj) {
591 bid = bid + contr1[nd + jj * contr1_dim1] *
592 hermit[ii + ((jj << 1) - 1) * hermit_dim1] +
593 contr2[nd + jj * contr2_dim1] * hermit[ii + (jj << 1) * hermit_dim1];
596 crvjac[ii + nd * crvjac_dim1] = bid;
602 /* ***********************************************************************
604 /* RETURN CALLING PROGRAM */
605 /* ***********************************************************************
609 AdvApp2Var_SysBase::mgsomsg_("MMA1CNT", 7L);
615 //=======================================================================
616 //function : mma1fdi_
618 //=======================================================================
619 int mma1fdi_(integer *ndimen,
621 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
635 /* System generated locals */
636 integer fpntab_dim1, somtab_dim1, somtab_offset, diftab_dim1,
637 diftab_offset, contr1_dim1, contr1_offset, contr2_dim1,
638 contr2_offset, i__1, i__2;
641 /* Local variables */
642 integer ideb, ifin, nroo2, ideru, iderv;
644 integer ii, nd, ibb, iim, nbp, iip;
645 doublereal bid1, bid2;
647 /* **********************************************************************
652 /* DiscretiZation of a non-polynomial function F(U,V) or of */
653 /* its derivative with fixed isoparameter. */
657 /* ALL, AB_SPECIFI::FONCTION&, DISCRETISATION, &POINT */
659 /* INPUT ARGUMENTS : */
660 /* ------------------ */
661 /* NDIMEN: Space dimension. */
662 /* UVFONC: Limits of the path of definition by U and by V of the approximated function */
663 /* FONCNP: The NAME of the non-polynomial function to be approximated */
664 /* (external program). */
665 /* ISOFAV: Fixed isoparameter for the discretization; */
666 /* = 1, discretization with fixed U and variable V. */
667 /* = 2, discretization with fixed V and variable U. */
668 /* TCONST: Iso value is also fixed. */
669 /* NBROOT: Number of INTERNAL discretization parameters. */
670 /* (if there are constraints, 2 extremities should be added).
672 /* This is also the root number of the Legendre polynom where */
673 /* the discretization is done. */
674 /* TTABLE: Table of discretization parameters and of 2 extremities */
675 /* (Respectively (-1, NBROOT Legendre roots,1) */
676 /* reframed within the adequate interval. */
677 /* IORDRE: Order of constraint imposed on the extremities of the iso. */
678 /* (If Iso-U, it is necessary to calculate the derivatives by V and vice */
680 /* = 0, the extremities of the iso are calculated. */
681 /* = 1, additionally the 1st derivative in the direction of the iso is calculated */
682 /* = 2, additionally the 2nd derivative in the direction of the iso is calculated */
683 /* IDERIV: Order of derivative transversal to fixed iso (If Iso-U=Uc */
684 /* is fixed, the derivative of order IDERIV is discretized by U of */
685 /* F(Uc,v). Same if iso-V is fixed). */
686 /* Varies from 0 (positioning) to 2 (2nd derivative). */
688 /* OUTPUT ARGUMENTS : */
689 /* ------------------- */
690 /* FPNTAB: Auxiliary table.
691 SOMTAB: Table of NBROOT/2 sums of 2 index points */
692 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
693 /* DIFTAB: Table of NBROOT/2 differences of 2 index points */
694 /* NBROOT-II+1 and II, for II = 1, NBROOT/2 */
695 /* CONTR1: Contains, if IORDRE>=0, values IORDRE+1 in TTABLE(0)
697 /* (1st extremity) of derivatives of F(Uc,Ve) or F(Ue,Vc), */
699 /* CONTR2: Contains, if IORDRE>=0, values IORDRE+1 in */
700 /* TTABLE(NBROOT+1) (2nd extremity) of: */
701 /* If ISOFAV=1, derived of order IDERIV by U, derived */
702 /* ordre 0 to IORDRE by V of F(Uc,Ve) or Uc=TCONST */
703 /* (fixed iso value) and Ve is the fixed extremity. */
704 /* If ISOFAV=2, derivative of order IDERIV by V, derivative */
705 /* of order 0 to IORDRE by U of F(Ue,Vc) or Vc=TCONST */
706 /* (fixed iso value) and Ue is the fixed extremity. */
707 /* IERCOD: Error code > 100; Pb in evaluation of FONCNP, */
708 /* the returned error code is equal to error code of FONCNP + 100. */
711 /* ---------------- */
713 /* REFERENCES CALLED : */
714 /* ----------------------- */
716 /* DESCRIPTION/NOTES/LIMITATIONS : */
717 /* ----------------------------------- */
718 /* The results of discretization are arranged in 2 tables */
719 /* SOMTAB and DIFTAB to earn time during the */
720 /* calculation of coefficients of the approximation curve. */
722 /* If NBROOT is uneven in SOMTAB(0,*) and DIFTAB(0,*) one stores */
723 /* the values of the median root of Legendre (0.D0 in (-1,1)). */
725 /* Function F(u,v) defined in UVFONC is reparameterized in */
726 /* (-1,1)x(-1,1). Then 1st and 2nd derivatives are renormalized. */
729 /* **********************************************************************
732 /* Name of the routine */
735 /* Parameter adjustments */
737 diftab_dim1 = *nbroot / 2 + 1;
738 diftab_offset = diftab_dim1;
739 diftab -= diftab_offset;
740 somtab_dim1 = *nbroot / 2 + 1;
741 somtab_offset = somtab_dim1;
742 somtab -= somtab_offset;
743 fpntab_dim1 = *ndimen;
745 contr2_dim1 = *ndimen;
746 contr2_offset = contr2_dim1 + 1;
747 contr2 -= contr2_offset;
748 contr1_dim1 = *ndimen;
749 contr1_offset = contr1_dim1 + 1;
750 contr1 -= contr1_offset;
753 ibb = AdvApp2Var_SysBase::mnfndeb_();
755 AdvApp2Var_SysBase::mgenmsg_("MMA1FDI", 7L);
759 /* --------------- Definition of the nb of points to calculate --------------
761 /* --> If constraints, the limits are also taken */
765 /* --> Otherwise, only Legendre roots (reframed) are used
771 /* --> Nb of point to calculate. */
772 nbp = ifin - ideb + 1;
775 /* --------------- Determination of the order of global derivation --------
777 /* --> ISOFAV takes only values 1 or 2. */
778 /* if Iso-U, derive by U of order IDERIV */
782 d__1 = (uvfonc[4] - uvfonc[3]) / 2.;
783 renor = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
784 /* if Iso-V, derive by V of order IDERIV */
788 d__1 = (uvfonc[6] - uvfonc[5]) / 2.;
789 renor = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
792 /* ----------- Discretization on roots of the ---------------
794 /* ---------------------- Legendre polynom of degree NBROOT -------------------
797 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (ndimen,
806 &fpntab[ideb * fpntab_dim1 + 1],
812 for (nd = 1; nd <= i__1; ++nd) {
814 for (ii = 1; ii <= i__2; ++ii) {
815 iip = (*nbroot + 1) / 2 + ii;
816 iim = nroo2 - ii + 1;
817 bid1 = fpntab[nd + iim * fpntab_dim1];
818 bid2 = fpntab[nd + iip * fpntab_dim1];
819 somtab[ii + nd * somtab_dim1] = renor * (bid2 + bid1);
820 diftab[ii + nd * diftab_dim1] = renor * (bid2 - bid1);
826 /* ------------ Case when discretisation is done on roots of a ----
828 /* ---------- Legendre polynom of uneven degree, 0 is root --------
831 if (*nbroot % 2 == 1) {
833 for (nd = 1; nd <= i__1; ++nd) {
834 somtab[nd * somtab_dim1] = renor * fpntab[nd + (nroo2 + 1) *
836 diftab[nd * diftab_dim1] = renor * fpntab[nd + (nroo2 + 1) *
842 for (nd = 1; nd <= i__1; ++nd) {
843 somtab[nd * somtab_dim1] = 0.;
844 diftab[nd * diftab_dim1] = 0.;
849 /* --------------------- Take into account constraints ----------------
853 /* --> Recover already calculated extremities. */
855 for (nd = 1; nd <= i__1; ++nd) {
856 contr1[nd + contr1_dim1] = renor * fpntab[nd];
857 contr2[nd + contr2_dim1] = renor * fpntab[nd + (*nbroot + 1) *
861 /* --> Nb of points to calculate/call to FONCNP */
863 /* If Iso-U, derive by V till order IORDRE */
865 /* --> Factor of normalisation 1st derivative. */
866 bid1 = (uvfonc[6] - uvfonc[5]) / 2.;
868 for (iderv = 1; iderv <= i__1; ++iderv) {
869 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
870 ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
871 nbp, ttable, &ideru, &iderv, &contr1[(iderv + 1) *
872 contr1_dim1 + 1], iercod);
879 for (iderv = 1; iderv <= i__1; ++iderv) {
880 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
881 ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
882 nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
883 iderv + 1) * contr2_dim1 + 1], iercod);
889 /* If Iso-V, derive by U till order IORDRE */
891 /* --> Factor of normalization 1st derivative. */
892 bid1 = (uvfonc[4] - uvfonc[3]) / 2.;
894 for (ideru = 1; ideru <= i__1; ++ideru) {
895 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
896 ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
897 nbp, ttable, &ideru, &iderv, &contr1[(ideru + 1) *
898 contr1_dim1 + 1], iercod);
905 for (ideru = 1; ideru <= i__1; ++ideru) {
906 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
907 ndimen, &uvfonc[3], &uvfonc[5], isofav, tconst, &
908 nbp, &ttable[*nbroot + 1], &ideru, &iderv, &contr2[(
909 ideru + 1) * contr2_dim1 + 1], iercod);
917 /* ------------------------- Normalization of derivatives -------------
919 /* (The function is redefined on (-1,1)*(-1,1)) */
922 for (ii = 1; ii <= i__1; ++ii) {
925 for (nd = 1; nd <= i__2; ++nd) {
926 contr1[nd + (ii + 1) * contr1_dim1] *= bid2;
927 contr2[nd + (ii + 1) * contr2_dim1] *= bid2;
934 /* ------------------------------ The end -------------------------------
940 AdvApp2Var_SysBase::maermsg_("MMA1FDI", iercod, 7L);
943 AdvApp2Var_SysBase::mgsomsg_("MMA1FDI", 7L);
948 //=======================================================================
949 //function : mma1fer_
951 //=======================================================================
952 int mma1fer_(integer *,//ndimen,
966 /* System generated locals */
967 integer crvjac_dim1, crvjac_offset, i__1, i__2;
969 /* Local variables */
970 integer idim, ncfja, ncfnw, ndses, ii, kk, ibb, ier;
974 /* ***********************************************************************
979 /* Calculate the degree and the errors of approximation of a border. */
983 /* TOUS,AB_SPECIFI :: COURBE&,TRONCATURE, &PRECISION */
985 /* INPUT ARGUMENTS : */
986 /* -------------------- */
988 /* NDIMEN: Total Dimension of the space (sum of dimensions of sub-spaces) */
989 /* NBSESP: Number of "independent" sub-spaces. */
990 /* NDIMSE: Table of dimensions of sub-spaces. */
991 /* IORDRE: Order of constraint at the extremities of the border */
992 /* -1 = no constraints, */
993 /* 0 = constraints of passage to limits (i.e. C0), */
994 /* 1 = C0 + constraintes of 1st derivatives (i.e. C1), */
995 /* 2 = C1 + constraintes of 2nd derivatives (i.e. C2). */
996 /* NDGJAC: Degree of development in series to use for the calculation */
997 /* in the base of Jacobi. */
998 /* CRVJAC: Table of coeff. of the curve of approximation in the */
999 /* base of Jacobi. */
1000 /* NCFLIM: Max number of coeff of the polynomial curve */
1001 /* of approximation (should be above or equal to */
1002 /* 2*IORDRE+2 and below or equal to 50). */
1003 /* EPSAPR: Table of errors of approximations that cannot be passed, */
1004 /* sub-space by sub-space. */
1006 /* OUTPUT ARGUMENTS : */
1007 /* --------------------- */
1008 /* YCVMAX: Auxiliary Table. */
1009 /* ERRMAX: Table of errors (sub-space by sub-space) */
1010 /* MAXIMUM made in the approximation of FONCNP by */
1012 /* ERRMOY: Table of errors (sub-space by sub-space) */
1013 /* AVERAGE made in the approximation of FONCNP by */
1015 /* NCOEFF: Number of significative coeffs. of the calculated "curve". */
1016 /* IERCOD: Error code */
1018 /* =-1, warning, required tolerance can't be */
1019 /* met with coefficients NFCLIM. */
1020 /* = 1, order of constraints (IORDRE) is not within authorised values */
1023 /* COMMONS USED : */
1024 /* ------------------ */
1026 /* REFERENCES CALLED : */
1027 /* --------------------- */
1029 /* DESCRIPTION/NOTES/LIMITATIONS : */
1030 /* ----------------------------------- */
1032 /* **********************************************************************
1035 /* Name of the routine */
1038 /* Parameter adjustments */
1044 crvjac_dim1 = *ndgjac + 1;
1045 crvjac_offset = crvjac_dim1;
1046 crvjac -= crvjac_offset;
1049 ibb = AdvApp2Var_SysBase::mnfndeb_();
1051 AdvApp2Var_SysBase::mgenmsg_("MMA1FER", 7L);
1056 ncfja = *ndgjac + 1;
1058 /* ------------ Calculate the degree of the curve and of the Max error --------
1060 /* -------------- of approximation for all sub-spaces --------
1064 for (ii = 1; ii <= i__1; ++ii) {
1067 /* ------------ cutting of coeff. and calculation of Max error -------
1070 AdvApp2Var_MathBase::mmtrpjj_(&ncfja, &ndses, &ncfja, &epsapr[ii], iordre, &crvjac[idim *
1071 crvjac_dim1], &ycvmax[1], &errmax[ii], &ncfnw);
1073 /* ******************************************************************
1075 /* ------------- If precision OK, calculate the average error -------
1077 /* ******************************************************************
1080 if (ncfnw <= *ncflim) {
1081 mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1082 crvjac_dim1], &ncfnw, &errmoy[ii]);
1083 *ncoeff = advapp_max(ncfnw,*ncoeff);
1085 /* ------------- Set the declined coefficients to 0.D0 -----------
1088 nbr0 = *ncflim - ncfnw;
1091 for (kk = 1; kk <= i__2; ++kk) {
1092 AdvApp2Var_SysBase::mvriraz_(&nbr0,
1093 &crvjac[ncfnw + (idim + kk - 1) * crvjac_dim1]);
1099 /* **************************************************************
1101 /* ------------------- If required precision can't be reached----
1103 /* **************************************************************
1108 /* ------------------------- calculate the Max error ------------
1111 AdvApp2Var_MathBase::mmaperx_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1112 crvjac_dim1], ncflim, &ycvmax[1], &errmax[ii], &ier);
1117 /* -------------------- nb of coeff to be returned -------------
1122 /* ------------------- and calculation of the average error ----
1125 mmaperm_(&ncfja, &ndses, &ncfja, iordre, &crvjac[idim *
1126 crvjac_dim1], ncflim, &errmoy[ii]);
1134 /* ------------------------------ The end -------------------------------
1136 /* --> The order of constraints is not within autorized values. */
1143 AdvApp2Var_SysBase::maermsg_("MMA1FER", iercod, 7L);
1146 AdvApp2Var_SysBase::mgsomsg_("MMA1FER", 7L);
1152 //=======================================================================
1153 //function : mma1her_
1155 //=======================================================================
1156 int AdvApp2Var_ApproxF2var::mma1her_(const integer *iordre,
1160 /* System generated locals */
1161 integer hermit_dim1, hermit_offset;
1163 /* Local variables */
1168 /* **********************************************************************
1173 /* Calculate 2*(IORDRE+1) Hermit polynoms of degree 2*IORDRE+1 */
1178 /* ALL, AB_SPECIFI::CONTRAINTE&, INTERPOLATION, &POLYNOME */
1180 /* INPUT ARGUMENTS : */
1181 /* ------------------ */
1182 /* IORDRE: Order of constraint. */
1183 /* = 0, Polynom of interpolation of order C0 on (-1,1). */
1184 /* = 1, Polynom of interpolation of order C0 and C1 on (-1,1). */
1185 /* = 2, Polynom of interpolation of order C0, C1 and C2 on (-1,1).
1188 /* OUTPUT ARGUMENTS : */
1189 /* ------------------- */
1190 /* HERMIT: Table of 2*IORDRE+2 coeff. of each of 2*(IORDRE+1) */
1191 /* HERMIT polynom. */
1192 /* IERCOD: Error code, */
1194 /* = 1, required order of constraint is not managed here. */
1195 /* COMMONS USED : */
1196 /* ---------------- */
1198 /* REFERENCES CALLED : */
1199 /* ----------------------- */
1201 /* DESCRIPTION/NOTES/LIMITATIONS : */
1202 /* ----------------------------------- */
1203 /* The part of HERMIT(*,2*i+j) table where j=1 or 2 and i=0 to IORDRE, */
1204 /* contains the coefficients of the polynom of degree 2*IORDRE+1 */
1205 /* such as ALL values in -1 and in +1 of this polynom and its */
1206 /* derivatives till order of derivation IORDRE are NULL, */
1207 /* EXCEPT for the derivative of order i: */
1208 /* - valued 1 in -1 if j=1 */
1209 /* - valued 1 in +1 if j=2. */
1211 /* **********************************************************************
1214 /* Name of the routine */
1217 /* Parameter adjustments */
1218 hermit_dim1 = (*iordre + 1) << 1;
1219 hermit_offset = hermit_dim1 + 1;
1220 hermit -= hermit_offset;
1223 ibb = AdvApp2Var_SysBase::mnfndeb_();
1225 AdvApp2Var_SysBase::mgenmsg_("MMA1HER", 7L);
1229 /* --- Recover (IORDRE+2) coeff of 2*(IORDRE+1) Hermit polynoms --
1233 hermit[hermit_dim1 + 1] = .5;
1234 hermit[hermit_dim1 + 2] = -.5;
1236 hermit[(hermit_dim1 << 1) + 1] = .5;
1237 hermit[(hermit_dim1 << 1) + 2] = .5;
1238 } else if (*iordre == 1) {
1239 hermit[hermit_dim1 + 1] = .5;
1240 hermit[hermit_dim1 + 2] = -.75;
1241 hermit[hermit_dim1 + 3] = 0.;
1242 hermit[hermit_dim1 + 4] = .25;
1244 hermit[(hermit_dim1 << 1) + 1] = .5;
1245 hermit[(hermit_dim1 << 1) + 2] = .75;
1246 hermit[(hermit_dim1 << 1) + 3] = 0.;
1247 hermit[(hermit_dim1 << 1) + 4] = -.25;
1249 hermit[hermit_dim1 * 3 + 1] = .25;
1250 hermit[hermit_dim1 * 3 + 2] = -.25;
1251 hermit[hermit_dim1 * 3 + 3] = -.25;
1252 hermit[hermit_dim1 * 3 + 4] = .25;
1254 hermit[(hermit_dim1 << 2) + 1] = -.25;
1255 hermit[(hermit_dim1 << 2) + 2] = -.25;
1256 hermit[(hermit_dim1 << 2) + 3] = .25;
1257 hermit[(hermit_dim1 << 2) + 4] = .25;
1258 } else if (*iordre == 2) {
1259 hermit[hermit_dim1 + 1] = .5;
1260 hermit[hermit_dim1 + 2] = -.9375;
1261 hermit[hermit_dim1 + 3] = 0.;
1262 hermit[hermit_dim1 + 4] = .625;
1263 hermit[hermit_dim1 + 5] = 0.;
1264 hermit[hermit_dim1 + 6] = -.1875;
1266 hermit[(hermit_dim1 << 1) + 1] = .5;
1267 hermit[(hermit_dim1 << 1) + 2] = .9375;
1268 hermit[(hermit_dim1 << 1) + 3] = 0.;
1269 hermit[(hermit_dim1 << 1) + 4] = -.625;
1270 hermit[(hermit_dim1 << 1) + 5] = 0.;
1271 hermit[(hermit_dim1 << 1) + 6] = .1875;
1273 hermit[hermit_dim1 * 3 + 1] = .3125;
1274 hermit[hermit_dim1 * 3 + 2] = -.4375;
1275 hermit[hermit_dim1 * 3 + 3] = -.375;
1276 hermit[hermit_dim1 * 3 + 4] = .625;
1277 hermit[hermit_dim1 * 3 + 5] = .0625;
1278 hermit[hermit_dim1 * 3 + 6] = -.1875;
1280 hermit[(hermit_dim1 << 2) + 1] = -.3125;
1281 hermit[(hermit_dim1 << 2) + 2] = -.4375;
1282 hermit[(hermit_dim1 << 2) + 3] = .375;
1283 hermit[(hermit_dim1 << 2) + 4] = .625;
1284 hermit[(hermit_dim1 << 2) + 5] = -.0625;
1285 hermit[(hermit_dim1 << 2) + 6] = -.1875;
1287 hermit[hermit_dim1 * 5 + 1] = .0625;
1288 hermit[hermit_dim1 * 5 + 2] = -.0625;
1289 hermit[hermit_dim1 * 5 + 3] = -.125;
1290 hermit[hermit_dim1 * 5 + 4] = .125;
1291 hermit[hermit_dim1 * 5 + 5] = .0625;
1292 hermit[hermit_dim1 * 5 + 6] = -.0625;
1294 hermit[hermit_dim1 * 6 + 1] = .0625;
1295 hermit[hermit_dim1 * 6 + 2] = .0625;
1296 hermit[hermit_dim1 * 6 + 3] = -.125;
1297 hermit[hermit_dim1 * 6 + 4] = -.125;
1298 hermit[hermit_dim1 * 6 + 5] = .0625;
1299 hermit[hermit_dim1 * 6 + 6] = .0625;
1304 /* ------------------------------ The End -------------------------------
1307 AdvApp2Var_SysBase::maermsg_("MMA1HER", iercod, 7L);
1309 AdvApp2Var_SysBase::mgsomsg_("MMA1HER", 7L);
1313 //=======================================================================
1314 //function : mma1jak_
1316 //=======================================================================
1317 int mma1jak_(integer *ndimen,
1327 /* System generated locals */
1328 integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
1329 crvjac_dim1, crvjac_offset, cgauss_dim1;
1331 /* Local variables */
1334 /* **********************************************************************
1339 /* Calculate the curve of approximation of a non-polynomial function */
1340 /* in the base of Jacobi. */
1344 /* FUNCTION,DISCRETISATION,APPROXIMATION,CONSTRAINT,CURVE,JACOBI */
1346 /* INPUT ARGUMENTS : */
1347 /* ------------------ */
1348 /* NDIMEN: Total dimension of the space (sum of dimensions */
1349 /* of sub-spaces) */
1350 /* NBROOT: Nb of points of discretization of the iso, extremities not */
1352 /* IORDRE: Order of constraint at the extremities of the boundary */
1353 /* -1 = no constraints, */
1354 /* 0 = constraints of passage of limits (i.e. C0), */
1355 /* 1 = C0 + constraints of 1st derivatives (i.e. C1), */
1356 /* 2 = C1 + constraints of 2nd derivatives (i.e. C2). */
1357 /* NDGJAC: Degree of development in series to be used for calculation in the */
1358 /* base of Jacobi. */
1360 /* OUTPUT ARGUMENTS : */
1361 /* ------------------- */
1362 /* CRVJAC : Curve of approximation of FONCNP with (eventually) */
1363 /* taking into account of constraints at the extremities. */
1364 /* This curve is of degree NDGJAC. */
1365 /* IERCOD : Error code : */
1366 /* 0 = All is ok. */
1367 /* 33 = Pb to return data of du block data */
1368 /* of coeff. of integration by GAUSS method. */
1369 /* by program MMAPPTT. */
1371 /* COMMONS USED : */
1372 /* ---------------- */
1374 /* REFERENCES CALLED : */
1375 /* ----------------------- */
1377 /* **********************************************************************
1380 /* Name of the routine */
1382 /* Parameter adjustments */
1383 diftab_dim1 = *nbroot / 2 + 1;
1384 diftab_offset = diftab_dim1;
1385 diftab -= diftab_offset;
1386 somtab_dim1 = *nbroot / 2 + 1;
1387 somtab_offset = somtab_dim1;
1388 somtab -= somtab_offset;
1389 crvjac_dim1 = *ndgjac + 1;
1390 crvjac_offset = crvjac_dim1;
1391 crvjac -= crvjac_offset;
1392 cgauss_dim1 = *nbroot / 2 + 1;
1395 ibb = AdvApp2Var_SysBase::mnfndeb_();
1397 AdvApp2Var_SysBase::mgenmsg_("MMA1JAK", 7L);
1401 /* ----------------- Recover coeffs of integration by Gauss -----------
1404 AdvApp2Var_ApproxF2var::mmapptt_(ndgjac, nbroot, iordre, cgauss, iercod);
1410 /* --------------- Calculate the curve in the base of Jacobi -----------
1413 mmmapcoe_(ndimen, ndgjac, iordre, nbroot, &somtab[somtab_offset], &diftab[
1414 diftab_offset], cgauss, &crvjac[crvjac_offset]);
1416 /* ------------------------------ The End -------------------------------
1421 AdvApp2Var_SysBase::maermsg_("MMA1JAK", iercod, 7L);
1424 AdvApp2Var_SysBase::mgsomsg_("MMA1JAK", 7L);
1429 //=======================================================================
1430 //function : mma1noc_
1432 //=======================================================================
1433 int mma1noc_(doublereal *dfuvin,
1442 /* System generated locals */
1446 /* Local variables */
1447 doublereal rider, riord;
1450 /* **********************************************************************
1455 /* Normalization of constraints of derivatives, defined on DFUVIN */
1456 /* on block DUVOUT. */
1460 /* ALL, AB_SPECIFI::VECTEUR&,DERIVEE&,NORMALISATION,&VECTEUR */
1462 /* INPUT ARGUMENTS : */
1463 /* ------------------ */
1464 /* DFUVIN: Limits of the block of definition by U and by V where
1466 /* constraints CNTRIN are defined. */
1467 /* NDIMEN: Dimension of the space. */
1468 /* IORDRE: Order of constraint imposed at the extremities of the iso. */
1469 /* (if Iso-U, it is necessary to calculate derivatives by V and vice */
1471 /* = 0, the extremities of the iso are calculated */
1472 /* = 1, additionally the 1st derivative in the direction */
1473 /* of the iso is calculated */
1474 /* = 2, additionally the 2nd derivative in the direction */
1475 /* of the iso is calculated */
1476 /* CNTRIN: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1477 /* of order IORDRE of F(Uc,v) or of F(u,Vc), following the */
1478 /* value of ISOFAV, RENORMALIZED by u and v in (-1,1). */
1479 /* DUVOUT: Limits of the block of definition by U and by V where the */
1480 /* constraints CNTOUT will be defined. */
1481 /* ISOFAV: Isoparameter fixed for the discretization; */
1482 /* = 1, discretization with fixed U=Uc and variable V. */
1483 /* = 2, discretization with fixed V=Vc and variable U. */
1484 /* IDERIV: Ordre de derivee transverse a l'iso fixee (Si Iso-U=Uc */
1485 /* is fixed, the derivative of order IDERIV is discretized by U */
1486 /* of F(Uc,v). The same if iso-V is fixed). */
1487 /* Varies from (positioning) to 2 (2nd derivative). */
1489 /* OUTPUT ARGUMENTS : */
1490 /* ------------------- */
1491 /* CNTOUT: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1492 /* of order IORDRE of F(Uc,v) or of F(u,Vc), depending on the */
1493 /* value of ISOFAV, RENORMALIZED for u and v in DUVOUT. */
1495 /* COMMONS USED : */
1496 /* ---------------- */
1498 /* REFERENCES CALLED : */
1499 /* --------------------- */
1501 /* DESCRIPTION/NOTES/LIMITATIONS : */
1502 /* ------------------------------- */
1503 /* CNTRIN can be an output/input argument, */
1506 /* CALL MMA1NOC(DFUVIN,NDIMEN,IORDRE,CNTRIN,DUVOUT */
1507 /* 1 ,ISOFAV,IDERIV,CNTRIN) */
1511 /* **********************************************************************
1514 /* Name of the routine */
1517 /* Parameter adjustments */
1524 ibb = AdvApp2Var_SysBase::mnfndeb_();
1526 AdvApp2Var_SysBase::mgenmsg_("MMA1NOC", 7L);
1529 /* --------------- Determination of coefficients of normalization -------
1533 d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1534 rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1535 d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1536 riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1539 d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1540 rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1541 d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1542 riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1545 /* ------------- Renormalization of the vector of constraint ---------------
1548 bid = rider * riord;
1550 for (nd = 1; nd <= i__1; ++nd) {
1551 cntout[nd] = bid * cntrin[nd];
1555 /* ------------------------------ The end -------------------------------
1559 AdvApp2Var_SysBase::mgsomsg_("MMA1NOC", 7L);
1564 //=======================================================================
1565 //function : mma1nop_
1567 //=======================================================================
1568 int mma1nop_(integer *nbroot,
1576 /* System generated locals */
1579 /* Local variables */
1580 doublereal alinu, blinu, alinv, blinv;
1583 /* ***********************************************************************
1588 /* Normalization of parameters of an iso, starting from */
1589 /* parametric block and parameters on (-1,1). */
1593 /* TOUS,AB_SPECIFI :: ISO&,POINT&,NORMALISATION,&POINT,&ISO */
1595 /* INPUT ARGUMENTS : */
1596 /* -------------------- */
1597 /* NBROOT: Nb of points of discretisation INSIDE the iso */
1598 /* defined on (-1,1). */
1599 /* ROOTLG: Table of discretization parameters on )-1,1( */
1601 /* UVFONC: Block of definition of the iso */
1602 /* ISOFAV: = 1, this is iso-u; =2, this is iso-v. */
1604 /* OUTPUT ARGUMENTS : */
1605 /* --------------------- */
1606 /* TTABLE: Table of parameters renormalized on UVFONC of the iso.
1608 /* IERCOD: = 0, OK */
1609 /* = 1, ISOFAV is out of allowed values. */
1612 /* **********************************************************************
1614 /* Name of the routine */
1617 /* Parameter adjustments */
1622 ibb = AdvApp2Var_SysBase::mnfndeb_();
1624 AdvApp2Var_SysBase::mgenmsg_("MMA1NOP", 7L);
1627 alinu = (uvfonc[4] - uvfonc[3]) / 2.;
1628 blinu = (uvfonc[4] + uvfonc[3]) / 2.;
1629 alinv = (uvfonc[6] - uvfonc[5]) / 2.;
1630 blinv = (uvfonc[6] + uvfonc[5]) / 2.;
1633 ttable[0] = uvfonc[5];
1635 for (ii = 1; ii <= i__1; ++ii) {
1636 ttable[ii] = alinv * rootlg[ii] + blinv;
1639 ttable[*nbroot + 1] = uvfonc[6];
1640 } else if (*isofav == 2) {
1641 ttable[0] = uvfonc[3];
1643 for (ii = 1; ii <= i__1; ++ii) {
1644 ttable[ii] = alinu * rootlg[ii] + blinu;
1647 ttable[*nbroot + 1] = uvfonc[4];
1654 /* ------------------------------ THE END -------------------------------
1663 AdvApp2Var_SysBase::maermsg_("MMA1NOP", iercod, 7L);
1666 AdvApp2Var_SysBase::mgsomsg_("MMA1NOP", 7L);
1673 //=======================================================================
1674 //function : mma2ac1_
1676 //=======================================================================
1677 int AdvApp2Var_ApproxF2var::mma2ac1_(integer const *ndimen,
1678 integer const *mxujac,
1679 integer const *mxvjac,
1680 integer const *iordru,
1681 integer const *iordrv,
1682 doublereal const *contr1,
1683 doublereal const * contr2,
1684 doublereal const *contr3,
1685 doublereal const *contr4,
1686 doublereal const *uhermt,
1687 doublereal const *vhermt,
1691 /* System generated locals */
1692 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
1693 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
1694 contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1,
1695 uhermt_offset, vhermt_dim1, vhermt_offset, patjac_dim1,
1696 patjac_dim2, patjac_offset, i__1, i__2, i__3, i__4, i__5;
1698 /* Local variables */
1701 doublereal bidu1, bidu2, bidv1, bidv2;
1702 integer ioru1, iorv1, ii, nd, jj, ku, kv;
1703 doublereal cnt1, cnt2, cnt3, cnt4;
1705 /* **********************************************************************
1710 /* Add polynoms of edge constraints. */
1714 /* TOUS,AB_SPECIFI::POINT&,CONTRAINTE&,ADDITION,&POLYNOME */
1716 /* INPUT ARGUMENTS : */
1717 /* ------------------ */
1718 /* NDIMEN: Dimension of the space. */
1719 /* MXUJAC: Max degree of the polynom of approximation by U. The */
1720 /* representation in the orthogonal base starts from degree */
1721 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1722 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1723 /* MXVJAC: Max degree of the polynom of approximation by V. The */
1724 /* representation in the orthogonal base starts from degree */
1725 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1726 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1727 /* IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
1728 /* to the step of constraints: C0, C1 or C2. */
1729 /* IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1730 /* to the step of constraints: C0, C1 or C2. */
1731 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
1732 /* extremities of F(U0,V0) and its derivatives. */
1733 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
1734 /* extremities of F(U1,V0) and its derivatives. */
1735 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
1736 /* extremities of F(U0,V1) and its derivatives. */
1737 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
1738 /* extremities of F(U1,V1) and its derivatives. */
1739 /* UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
1740 /* VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1741 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1742 /* of F(u,v) WITHOUT taking into account the constraints. */
1744 /* OUTPUT ARGUMENTS : */
1745 /* ------------------- */
1746 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1747 /* of F(u,v) WITH taking into account of constraints. */
1749 /* **********************************************************************
1751 /* Name of the routine */
1753 /* --------------------------- Initialization --------------------------
1756 /* Parameter adjustments */
1757 patjac_dim1 = *mxujac + 1;
1758 patjac_dim2 = *mxvjac + 1;
1759 patjac_offset = patjac_dim1 * patjac_dim2;
1760 patjac -= patjac_offset;
1761 uhermt_dim1 = (*iordru << 1) + 2;
1762 uhermt_offset = uhermt_dim1;
1763 uhermt -= uhermt_offset;
1764 vhermt_dim1 = (*iordrv << 1) + 2;
1765 vhermt_offset = vhermt_dim1;
1766 vhermt -= vhermt_offset;
1767 contr4_dim1 = *ndimen;
1768 contr4_dim2 = *iordru + 2;
1769 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
1770 contr4 -= contr4_offset;
1771 contr3_dim1 = *ndimen;
1772 contr3_dim2 = *iordru + 2;
1773 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
1774 contr3 -= contr3_offset;
1775 contr2_dim1 = *ndimen;
1776 contr2_dim2 = *iordru + 2;
1777 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
1778 contr2 -= contr2_offset;
1779 contr1_dim1 = *ndimen;
1780 contr1_dim2 = *iordru + 2;
1781 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
1782 contr1 -= contr1_offset;
1785 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1787 AdvApp2Var_SysBase::mgenmsg_("MMA2AC1", 7L);
1790 /* ------------ SUBTRACTION OF ANGULAR CONSTRAINTS -------------------
1793 ioru1 = *iordru + 1;
1794 iorv1 = *iordrv + 1;
1795 ndgu = (*iordru << 1) + 1;
1796 ndgv = (*iordrv << 1) + 1;
1799 for (jj = 1; jj <= i__1; ++jj) {
1801 for (ii = 1; ii <= i__2; ++ii) {
1803 for (nd = 1; nd <= i__3; ++nd) {
1804 cnt1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
1805 cnt2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
1806 cnt3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
1807 cnt4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
1809 for (kv = 0; kv <= i__4; ++kv) {
1810 bidv1 = vhermt[kv + ((jj << 1) - 1) * vhermt_dim1];
1811 bidv2 = vhermt[kv + (jj << 1) * vhermt_dim1];
1813 for (ku = 0; ku <= i__5; ++ku) {
1814 bidu1 = uhermt[ku + ((ii << 1) - 1) * uhermt_dim1];
1815 bidu2 = uhermt[ku + (ii << 1) * uhermt_dim1];
1816 patjac[ku + (kv + nd * patjac_dim2) * patjac_dim1] =
1817 patjac[ku + (kv + nd * patjac_dim2) *
1818 patjac_dim1] - bidu1 * bidv1 * cnt1 - bidu2 *
1819 bidv1 * cnt2 - bidu1 * bidv2 * cnt3 - bidu2 *
1832 /* ------------------------------ The end -------------------------------
1836 AdvApp2Var_SysBase::mgsomsg_("MMA2AC1", 7L);
1841 //=======================================================================
1842 //function : mma2ac2_
1844 //=======================================================================
1845 int AdvApp2Var_ApproxF2var::mma2ac2_(const integer *ndimen,
1846 const integer *mxujac,
1847 const integer *mxvjac,
1848 const integer *iordrv,
1849 const integer *nclimu,
1850 const integer *ncfiv1,
1851 const doublereal *crbiv1,
1852 const integer *ncfiv2,
1853 const doublereal *crbiv2,
1854 const doublereal *vhermt,
1858 /* System generated locals */
1859 integer crbiv1_dim1, crbiv1_dim2, crbiv1_offset, crbiv2_dim1, crbiv2_dim2,
1860 crbiv2_offset, patjac_dim1, patjac_dim2, patjac_offset,
1861 vhermt_dim1, vhermt_offset, i__1, i__2, i__3, i__4;
1863 /* Local variables */
1865 integer ndgv1, ndgv2, ii, jj, nd, kk;
1866 doublereal bid1, bid2;
1868 /* **********************************************************************
1873 /* Add polynoms of constraints */
1877 /* FUNCTION,APPROXIMATION,COEFFICIENT,POLYNOM */
1879 /* INPUT ARGUMENTS : */
1880 /* ------------------ */
1881 /* NDIMEN: Dimension of the space. */
1882 /* MXUJAC: Max degree of the polynom of approximation by U. The */
1883 /* representation in the orthogonal base starts from degree */
1884 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1885 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1886 /* MXVJAC: Max degree of the polynom of approximation by V. The */
1887 /* representation in the orthogonal base starts from degree */
1888 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1889 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1890 /* IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1891 /* to the step of constraints: C0, C1 or C2. */
1892 /* NCLIMU LIMIT nb of coeff by u of the solution P(u,v)
1893 * NCFIV1: Nb of Coeff. of curves stored in CRBIV1. */
1894 /* CRBIV1: Table of coeffs of the approximation of iso-V0 and its */
1895 /* derivatives till order IORDRV. */
1896 /* NCFIV2: Nb of Coeff. of curves stored in CRBIV2. */
1897 /* CRBIV2: Table of coeffs of approximation of iso-V1 and its */
1898 /* derivatives till order IORDRV. */
1899 /* VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1900 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1901 /* of F(u,v) WITHOUT taking into account the constraints. */
1903 /* OUTPUT ARGUMENTS : */
1904 /* ------------------- */
1905 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1906 /* of F(u,v) WITH taking into account of constraints. */
1911 /* **********************************************************************
1913 /* Name of the routine */
1915 /* --------------------------- Initialisations --------------------------
1918 /* Parameter adjustments */
1919 patjac_dim1 = *mxujac + 1;
1920 patjac_dim2 = *mxvjac + 1;
1921 patjac_offset = patjac_dim1 * patjac_dim2;
1922 patjac -= patjac_offset;
1923 vhermt_dim1 = (*iordrv << 1) + 2;
1924 vhermt_offset = vhermt_dim1;
1925 vhermt -= vhermt_offset;
1928 crbiv2_dim1 = *nclimu;
1929 crbiv2_dim2 = *ndimen;
1930 crbiv2_offset = crbiv2_dim1 * (crbiv2_dim2 + 1);
1931 crbiv2 -= crbiv2_offset;
1932 crbiv1_dim1 = *nclimu;
1933 crbiv1_dim2 = *ndimen;
1934 crbiv1_offset = crbiv1_dim1 * (crbiv1_dim2 + 1);
1935 crbiv1 -= crbiv1_offset;
1938 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1940 AdvApp2Var_SysBase::mgenmsg_("MMA2AC2", 7L);
1943 /* ------------ ADDING of coeff by u of curves, by v of Hermit --------
1947 for (ii = 1; ii <= i__1; ++ii) {
1948 ndgv1 = ncfiv1[ii] - 1;
1949 ndgv2 = ncfiv2[ii] - 1;
1951 for (nd = 1; nd <= i__2; ++nd) {
1952 i__3 = (*iordrv << 1) + 1;
1953 for (jj = 0; jj <= i__3; ++jj) {
1954 bid1 = vhermt[jj + ((ii << 1) - 1) * vhermt_dim1];
1956 for (kk = 0; kk <= i__4; ++kk) {
1957 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
1958 bid1 * crbiv1[kk + (nd + ii * crbiv1_dim2) *
1962 bid2 = vhermt[jj + (ii << 1) * vhermt_dim1];
1964 for (kk = 0; kk <= i__4; ++kk) {
1965 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
1966 bid2 * crbiv2[kk + (nd + ii * crbiv2_dim2) *
1977 /* ------------------------------ The end -------------------------------
1981 AdvApp2Var_SysBase::mgsomsg_("MMA2AC2", 7L);
1987 //=======================================================================
1988 //function : mma2ac3_
1990 //=======================================================================
1991 int AdvApp2Var_ApproxF2var::mma2ac3_(const integer *ndimen,
1992 const integer *mxujac,
1993 const integer *mxvjac,
1994 const integer *iordru,
1995 const integer *nclimv,
1996 const integer *ncfiu1,
1997 const doublereal * crbiu1,
1998 const integer *ncfiu2,
1999 const doublereal *crbiu2,
2000 const doublereal *uhermt,
2004 /* System generated locals */
2005 integer crbiu1_dim1, crbiu1_dim2, crbiu1_offset, crbiu2_dim1, crbiu2_dim2,
2006 crbiu2_offset, patjac_dim1, patjac_dim2, patjac_offset,
2007 uhermt_dim1, uhermt_offset, i__1, i__2, i__3, i__4;
2009 /* Local variables */
2011 integer ndgu1, ndgu2, ii, jj, nd, kk;
2012 doublereal bid1, bid2;
2014 /* **********************************************************************
2019 /* Ajout des polynomes de contraintes */
2023 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
2025 /* INPUT ARGUMENTS : */
2026 /* ------------------ */
2027 /* NDIMEN: Dimension of the space. */
2028 /* MXUJAC: Max degree of the polynom of approximation by U. The */
2029 /* representation in the orthogonal base starts from degree */
2030 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2031 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2032 /* MXVJAC: Max degree of the polynom of approximation by V. The */
2033 /* representation in the orthogonal base starts from degree */
2034 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2035 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2036 /* IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
2037 /* to the step of constraints: C0, C1 or C2. */
2038 /* NCLIMV LIMIT nb of coeff by v of the solution P(u,v)
2039 * NCFIU1: Nb of Coeff. of curves stored in CRBIU1. */
2040 /* CRBIU1: Table of coeffs of the approximation of iso-U0 and its */
2041 /* derivatives till order IORDRU. */
2042 /* NCFIU2: Nb of Coeff. of curves stored in CRBIU2. */
2043 /* CRBIU2: Table of coeffs of approximation of iso-U1 and its */
2044 /* derivatives till order IORDRU */
2045 /* UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
2046 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
2047 /* of F(u,v) WITHOUT taking into account the constraints. */
2049 /* OUTPUT ARGUMENTS : */
2050 /* ------------------- */
2051 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
2052 /* of F(u,v) WITH taking into account of constraints. */
2056 /* **********************************************************************
2058 /* The name of the routine */
2060 /* --------------------------- Initializations --------------------------
2063 /* Parameter adjustments */
2064 patjac_dim1 = *mxujac + 1;
2065 patjac_dim2 = *mxvjac + 1;
2066 patjac_offset = patjac_dim1 * patjac_dim2;
2067 patjac -= patjac_offset;
2068 uhermt_dim1 = (*iordru << 1) + 2;
2069 uhermt_offset = uhermt_dim1;
2070 uhermt -= uhermt_offset;
2073 crbiu2_dim1 = *nclimv;
2074 crbiu2_dim2 = *ndimen;
2075 crbiu2_offset = crbiu2_dim1 * (crbiu2_dim2 + 1);
2076 crbiu2 -= crbiu2_offset;
2077 crbiu1_dim1 = *nclimv;
2078 crbiu1_dim2 = *ndimen;
2079 crbiu1_offset = crbiu1_dim1 * (crbiu1_dim2 + 1);
2080 crbiu1 -= crbiu1_offset;
2083 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
2085 AdvApp2Var_SysBase::mgenmsg_("MMA2AC3", 7L);
2088 /* ------------ ADDING of coeff by u of curves, by v of Hermit --------
2092 for (ii = 1; ii <= i__1; ++ii) {
2093 ndgu1 = ncfiu1[ii] - 1;
2094 ndgu2 = ncfiu2[ii] - 1;
2096 for (nd = 1; nd <= i__2; ++nd) {
2098 for (jj = 0; jj <= i__3; ++jj) {
2099 bid1 = crbiu1[jj + (nd + ii * crbiu1_dim2) * crbiu1_dim1];
2100 i__4 = (*iordru << 1) + 1;
2101 for (kk = 0; kk <= i__4; ++kk) {
2102 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
2103 bid1 * uhermt[kk + ((ii << 1) - 1) * uhermt_dim1];
2109 for (jj = 0; jj <= i__3; ++jj) {
2110 bid2 = crbiu2[jj + (nd + ii * crbiu2_dim2) * crbiu2_dim1];
2111 i__4 = (*iordru << 1) + 1;
2112 for (kk = 0; kk <= i__4; ++kk) {
2113 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
2114 bid2 * uhermt[kk + (ii << 1) * uhermt_dim1];
2125 /* ------------------------------ The end -------------------------------
2129 AdvApp2Var_SysBase::mgsomsg_("MMA2AC3", 7L);
2134 //=======================================================================
2135 //function : mma2can_
2137 //=======================================================================
2138 int AdvApp2Var_ApproxF2var::mma2can_(const integer *ncfmxu,
2139 const integer *ncfmxv,
2140 const integer *ndimen,
2141 const integer *iordru,
2142 const integer *iordrv,
2143 const integer *ncoefu,
2144 const integer *ncoefv,
2145 const doublereal *patjac,
2151 /* System generated locals */
2152 integer patjac_dim1, patjac_dim2, patjac_offset, patcan_dim1, patcan_dim2,
2153 patcan_offset, i__1, i__2;
2155 /* Local variables */
2157 integer ilon1, ilon2, ii, nd;
2159 /* **********************************************************************
2164 /* Change of Jacobi base to canonical (-1,1) and writing in a greater */
2169 /* ALL,AB_SPECIFI,CARREAU&,CONVERSION,JACOBI,CANNONIQUE,&CARREAU */
2171 /* INPUT ARGUMENTS : */
2172 /* -------------------- */
2173 /* NCFMXU: Dimension by U of resulting table PATCAN */
2174 /* NCFMXV: Dimension by V of resulting table PATCAN */
2175 /* NDIMEN: Dimension of the workspace. */
2176 /* IORDRU: Order of constraint by U */
2177 /* IORDRV: Order of constraint by V. */
2178 /* NCOEFU: Nb of coeff by U of square PATJAC */
2179 /* NCOEFV: Nb of coeff by V of square PATJAC */
2180 /* PATJAC: Square in the base of Jacobi of order IORDRU by U and */
2183 /* OUTPUT ARGUMENTS : */
2184 /* --------------------- */
2185 /* PATAUX: Auxiliary Table. */
2186 /* PATCAN: Table of coefficients in the canonic base. */
2187 /* IERCOD: Error code. */
2188 /* = 0, everything goes well, and all things are equal. */
2189 /* = 1, the program refuses to process with incorrect input arguments */
2192 /* COMMONS USED : */
2193 /* ------------------ */
2195 /* REFERENCES CALLED : */
2196 /* --------------------- */
2198 /* DESCRIPTION/NOTES/LIMITATIONS : */
2199 /* ----------------------------------- */
2201 /* **********************************************************************
2205 /* Parameter adjustments */
2206 patcan_dim1 = *ncfmxu;
2207 patcan_dim2 = *ncfmxv;
2208 patcan_offset = patcan_dim1 * (patcan_dim2 + 1) + 1;
2209 patcan -= patcan_offset;
2211 patjac_dim1 = *ncoefu;
2212 patjac_dim2 = *ncoefv;
2213 patjac_offset = patjac_dim1 * (patjac_dim2 + 1) + 1;
2214 patjac -= patjac_offset;
2217 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
2219 AdvApp2Var_SysBase::mgenmsg_("MMA2CAN", 7L);
2223 if (*iordru < -1 || *iordru > 2) {
2226 if (*iordrv < -1 || *iordrv > 2) {
2229 if (*ncoefu > *ncfmxu || *ncoefv > *ncfmxv) {
2233 /* --> Pass to canonic base (-1,1) */
2234 mmjacpt_(ndimen, ncoefu, ncoefv, iordru, iordrv, &patjac[patjac_offset], &
2235 pataux[1], &patcan[patcan_offset]);
2237 /* --> Write all in a greater table */
2238 AdvApp2Var_MathBase::mmfmca8_(ncoefu,
2244 &patcan[patcan_offset],
2245 &patcan[patcan_offset]);
2247 /* --> Complete with zeros the resulting table. */
2248 ilon1 = *ncfmxu - *ncoefu;
2249 ilon2 = *ncfmxu * (*ncfmxv - *ncoefv);
2251 for (nd = 1; nd <= i__1; ++nd) {
2254 for (ii = 1; ii <= i__2; ++ii) {
2255 AdvApp2Var_SysBase::mvriraz_(&ilon1,
2256 &patcan[*ncoefu + 1 + (ii + nd * patcan_dim2) * patcan_dim1]);
2261 AdvApp2Var_SysBase::mvriraz_(&ilon2,
2262 &patcan[(*ncoefv + 1 + nd * patcan_dim2) * patcan_dim1 + 1]);
2269 /* ----------------------
2277 AdvApp2Var_SysBase::maermsg_("MMA2CAN", iercod, 7L);
2279 AdvApp2Var_SysBase::mgsomsg_("MMA2CAN", 7L);
2284 //=======================================================================
2285 //function : mma2cd1_
2287 //=======================================================================
2288 int mma2cd1_(integer *ndimen,
2311 /* System generated locals */
2312 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
2313 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
2314 contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1,
2315 uhermt_offset, vhermt_dim1, vhermt_offset, fpntbu_dim1,
2316 fpntbu_offset, fpntbv_dim1, fpntbv_offset, sosotb_dim1,
2317 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2318 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2319 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4,
2322 /* Local variables */
2323 integer ncfhu, ncfhv, nuroo, nvroo, nd, ii, jj, kk, ll, ibb, kkm,
2325 doublereal bid1, bid2, bid3, bid4;
2326 doublereal diu1, diu2, div1, div2, sou1, sou2, sov1, sov2;
2328 /* **********************************************************************
2333 /* Discretisation on the parameters of polynoms of interpolation */
2334 /* of constraints at the corners of order IORDRE. */
2338 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2340 /* INPUT ARGUMENTS : */
2341 /* ------------------ */
2342 /* NDIMEN: Dimension of the space. */
2343 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2344 /* This is also the nb of root of Legendre polynom where discretization is done. */
2345 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
2347 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
2348 /* This is also the nb of root of Legendre polynom where discretization is done. */
2349 /* VROOTL: Table of discretization parameters on (-1,1) by V. */
2350 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
2351 /* = 0, calculate the extremities of iso-V */
2352 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2353 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2354 /* IORDRV: Order of constraint imposed at the extremities of iso-U */
2355 /* = 0, calculate the extremities of iso-U */
2356 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
2357 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
2358 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
2359 /* extremities of F(U0,V0) and its derivatives. */
2360 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
2361 /* extremities of F(U1,V0) and its derivatives. */
2362 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
2363 /* extremities of F(U0,V1) and its derivatives. */
2364 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
2365 /* extremities of F(U1,V1) and its derivatives. */
2366 /* SOSOTB: Preinitialized table (input/output argument). */
2367 /* DISOTB: Preinitialized table (input/output argument). */
2368 /* SODITB: Preinitialized table (input/output argument). */
2369 /* DIDITB: Preinitialized table (input/output argument) */
2371 /* OUTPUT ARGUMENTS : */
2372 /* ------------------- */
2373 /* FPNTBU: Auxiliary table. */
2374 /* FPNTBV: Auxiliary table. */
2375 /* UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
2376 /* VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2377 /* SOSOTB: Table where the terms of constraints are added */
2378 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2379 /* with ui and vj positive roots of the Legendre polynom */
2380 /* of degree NBPNTU and NBPNTV respectively. */
2381 /* DISOTB: Table where the terms of constraints are added */
2382 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2383 /* with ui and vj positive roots of the polynom of Legendre */
2384 /* of degree NBPNTU and NBPNTV respectively. */
2385 /* SODITB: Table where the terms of constraints are added */
2386 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2387 /* with ui and vj positive roots of the polynom of Legendre */
2388 /* of degree NBPNTU and NBPNTV respectively. */
2389 /* DIDITB: Table where the terms of constraints are added */
2390 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2391 /* with ui and vj positive roots of the polynom of Legendre */
2392 /* of degree NBPNTU and NBPNTV respectively. */
2394 /* COMMONS USED : */
2395 /* ---------------- */
2397 /* REFERENCES CALLED : */
2398 /* ----------------------- */
2400 /* DESCRIPTION/NOTES/LIMITATIONS : */
2401 /* ----------------------------------- */
2404 /* **********************************************************************
2407 /* Name of the routine */
2410 /* Parameter adjustments */
2412 diditb_dim1 = *nbpntu / 2 + 1;
2413 diditb_dim2 = *nbpntv / 2 + 1;
2414 diditb_offset = diditb_dim1 * diditb_dim2;
2415 diditb -= diditb_offset;
2416 disotb_dim1 = *nbpntu / 2;
2417 disotb_dim2 = *nbpntv / 2;
2418 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2419 disotb -= disotb_offset;
2420 soditb_dim1 = *nbpntu / 2;
2421 soditb_dim2 = *nbpntv / 2;
2422 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2423 soditb -= soditb_offset;
2424 sosotb_dim1 = *nbpntu / 2 + 1;
2425 sosotb_dim2 = *nbpntv / 2 + 1;
2426 sosotb_offset = sosotb_dim1 * sosotb_dim2;
2427 sosotb -= sosotb_offset;
2429 uhermt_dim1 = (*iordru << 1) + 2;
2430 uhermt_offset = uhermt_dim1;
2431 uhermt -= uhermt_offset;
2432 fpntbu_dim1 = *nbpntu;
2433 fpntbu_offset = fpntbu_dim1 + 1;
2434 fpntbu -= fpntbu_offset;
2435 vhermt_dim1 = (*iordrv << 1) + 2;
2436 vhermt_offset = vhermt_dim1;
2437 vhermt -= vhermt_offset;
2438 fpntbv_dim1 = *nbpntv;
2439 fpntbv_offset = fpntbv_dim1 + 1;
2440 fpntbv -= fpntbv_offset;
2441 contr4_dim1 = *ndimen;
2442 contr4_dim2 = *iordru + 2;
2443 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
2444 contr4 -= contr4_offset;
2445 contr3_dim1 = *ndimen;
2446 contr3_dim2 = *iordru + 2;
2447 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
2448 contr3 -= contr3_offset;
2449 contr2_dim1 = *ndimen;
2450 contr2_dim2 = *iordru + 2;
2451 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
2452 contr2 -= contr2_offset;
2453 contr1_dim1 = *ndimen;
2454 contr1_dim2 = *iordru + 2;
2455 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
2456 contr1 -= contr1_offset;
2459 ibb = AdvApp2Var_SysBase::mnfndeb_();
2461 AdvApp2Var_SysBase::mgenmsg_("MMA2CD1", 7L);
2464 /* ------------------- Discretisation of Hermite polynoms -----------
2467 ncfhu = (*iordru + 1) << 1;
2469 for (ii = 1; ii <= i__1; ++ii) {
2471 for (ll = 1; ll <= i__2; ++ll) {
2472 AdvApp2Var_MathBase::mmmpocur_(&ncfhu, &c__1, &ncfhu, &uhermt[ii * uhermt_dim1], &
2473 urootl[ll], &fpntbu[ll + ii * fpntbu_dim1]);
2478 ncfhv = (*iordrv + 1) << 1;
2480 for (jj = 1; jj <= i__1; ++jj) {
2482 for (kk = 1; kk <= i__2; ++kk) {
2483 AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[jj * vhermt_dim1], &
2484 vrootl[kk], &fpntbv[kk + jj * fpntbv_dim1]);
2490 /* ---- The discretizations of polynoms of constraints are subtracted ----
2493 nuroo = *nbpntu / 2;
2494 nvroo = *nbpntv / 2;
2496 for (nd = 1; nd <= i__1; ++nd) {
2499 for (jj = 1; jj <= i__2; ++jj) {
2501 for (ii = 1; ii <= i__3; ++ii) {
2502 bid1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
2503 bid2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
2504 bid3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
2505 bid4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
2508 for (kk = 1; kk <= i__4; ++kk) {
2509 kkp = (*nbpntv + 1) / 2 + kk;
2510 kkm = nvroo - kk + 1;
2511 sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] +
2512 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2513 div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] -
2514 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2515 sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[kkm
2516 + (jj << 1) * fpntbv_dim1];
2517 div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[kkm
2518 + (jj << 1) * fpntbv_dim1];
2520 for (ll = 1; ll <= i__5; ++ll) {
2521 llp = (*nbpntu + 1) / 2 + ll;
2522 llm = nuroo - ll + 1;
2523 sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] +
2524 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2525 diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] -
2526 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2527 sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2528 llm + (ii << 1) * fpntbu_dim1];
2529 diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2530 llm + (ii << 1) * fpntbu_dim1];
2531 sosotb[ll + (kk + nd * sosotb_dim2) * sosotb_dim1] =
2532 sosotb[ll + (kk + nd * sosotb_dim2) *
2533 sosotb_dim1] - bid1 * sou1 * sov1 - bid2 *
2534 sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2536 soditb[ll + (kk + nd * soditb_dim2) * soditb_dim1] =
2537 soditb[ll + (kk + nd * soditb_dim2) *
2538 soditb_dim1] - bid1 * sou1 * div1 - bid2 *
2539 sou2 * div1 - bid3 * sou1 * div2 - bid4 *
2541 disotb[ll + (kk + nd * disotb_dim2) * disotb_dim1] =
2542 disotb[ll + (kk + nd * disotb_dim2) *
2543 disotb_dim1] - bid1 * diu1 * sov1 - bid2 *
2544 diu2 * sov1 - bid3 * diu1 * sov2 - bid4 *
2546 diditb[ll + (kk + nd * diditb_dim2) * diditb_dim1] =
2547 diditb[ll + (kk + nd * diditb_dim2) *
2548 diditb_dim1] - bid1 * diu1 * div1 - bid2 *
2549 diu2 * div1 - bid3 * diu1 * div2 - bid4 *
2556 /* ------------ Case when the discretization is done only on the roots
2558 /* ---------- of Legendre polynom of uneven degree, 0 is root
2561 if (*nbpntu % 2 == 1) {
2562 sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2563 sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2565 for (kk = 1; kk <= i__4; ++kk) {
2566 kkp = (*nbpntv + 1) / 2 + kk;
2567 kkm = nvroo - kk + 1;
2568 sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] +
2569 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2570 div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] -
2571 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2572 sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[
2573 kkm + (jj << 1) * fpntbv_dim1];
2574 div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[
2575 kkm + (jj << 1) * fpntbv_dim1];
2576 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1] =
2577 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1]
2578 - bid1 * sou1 * sov1 - bid2 * sou2 * sov1 -
2579 bid3 * sou1 * sov2 - bid4 * sou2 * sov2;
2580 diditb[(kk + nd * diditb_dim2) * diditb_dim1] =
2581 diditb[(kk + nd * diditb_dim2) * diditb_dim1]
2582 - bid1 * sou1 * div1 - bid2 * sou2 * div1 -
2583 bid3 * sou1 * div2 - bid4 * sou2 * div2;
2588 if (*nbpntv % 2 == 1) {
2589 sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2590 sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2592 for (ll = 1; ll <= i__4; ++ll) {
2593 llp = (*nbpntu + 1) / 2 + ll;
2594 llm = nuroo - ll + 1;
2595 sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] +
2596 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2597 diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] -
2598 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2599 sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2600 llm + (ii << 1) * fpntbu_dim1];
2601 diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2602 llm + (ii << 1) * fpntbu_dim1];
2603 sosotb[ll + nd * sosotb_dim2 * sosotb_dim1] = sosotb[
2604 ll + nd * sosotb_dim2 * sosotb_dim1] - bid1 *
2605 sou1 * sov1 - bid2 * sou2 * sov1 - bid3 *
2606 sou1 * sov2 - bid4 * sou2 * sov2;
2607 diditb[ll + nd * diditb_dim2 * diditb_dim1] = diditb[
2608 ll + nd * diditb_dim2 * diditb_dim1] - bid1 *
2609 diu1 * sov1 - bid2 * diu2 * sov1 - bid3 *
2610 diu1 * sov2 - bid4 * diu2 * sov2;
2615 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2616 sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2617 sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2618 sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2619 sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2620 sosotb[nd * sosotb_dim2 * sosotb_dim1] = sosotb[nd *
2621 sosotb_dim2 * sosotb_dim1] - bid1 * sou1 * sov1 -
2622 bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2624 diditb[nd * diditb_dim2 * diditb_dim1] = diditb[nd *
2625 diditb_dim2 * diditb_dim1] - bid1 * sou1 * sov1 -
2626 bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2638 /* ------------------------------ The End -------------------------------
2643 AdvApp2Var_SysBase::mgsomsg_("MMA2CD1", 7L);
2648 //=======================================================================
2649 //function : mma2cd2_
2651 //=======================================================================
2652 int mma2cd2_(integer *ndimen,
2670 /* System generated locals */
2671 integer sotbv1_dim1, sotbv1_dim2, sotbv1_offset, sotbv2_dim1, sotbv2_dim2,
2672 sotbv2_offset, ditbv1_dim1, ditbv1_dim2, ditbv1_offset,
2673 ditbv2_dim1, ditbv2_dim2, ditbv2_offset, fpntab_dim1,
2674 fpntab_offset, vhermt_dim1, vhermt_offset, sosotb_dim1,
2675 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2676 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2677 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2679 /* Local variables */
2680 integer ncfhv, nuroo, nvroo, ii, nd, jj, kk, ibb, jjm, jjp;
2681 doublereal bid1, bid2, bid3, bid4;
2683 /* **********************************************************************
2687 /* Discretisation on the parameters of polynoms of interpolation */
2688 /* of constraints on 2 borders iso-V of order IORDRV. */
2693 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2697 /* INPUT ARGUMENTS : */
2698 /* ------------------ */
2699 /* NDIMEN: Dimension of the space. */
2700 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2701 /* This is also the nb of root of Legendre polynom where discretization is done. */
2702 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
2704 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
2705 /* This is also the nb of root of Legendre polynom where discretization is done. */
2706 /* VROOTL: Table of discretization parameters on (-1,1) by V. */
2707 /* IORDRV: Order of constraint imposed at the extremities of iso-V */
2708 /* = 0, calculate the extremities of iso-V */
2709 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2710 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2711 /* SOTBV1: Table of NBPNTV/2 sums of 2 index points */
2712 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2713 /* SOTBV2: Table of NBPNTV/2 sums of 2 index points */
2714 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2715 /* DITBV1: Table of NBPNTV/2 differences of 2 index points */
2716 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2717 /* DITBV2: Table of NBPNTV/2 differences of 2 index points */
2718 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2719 /* SOSOTB: Preinitialized table (input/output argument). */
2720 /* DISOTB: Preinitialized table (input/output argument). */
2721 /* SODITB: Preinitialized table (input/output argument). */
2722 /* DIDITB: Preinitialized table (input/output argument) */
2724 /* OUTPUT ARGUMENTS : */
2725 /* ------------------- */
2726 /* FPNTAB: Auxiliary table. */
2727 /* VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2728 /* SOSOTB: Table where the terms of constraints are added */
2729 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2730 /* with ui and vj positive roots of the Legendre polynom */
2731 /* of degree NBPNTU and NBPNTV respectively. */
2732 /* DISOTB: Table where the terms of constraints are added */
2733 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2734 /* with ui and vj positive roots of the polynom of Legendre */
2735 /* of degree NBPNTU and NBPNTV respectively. */
2736 /* SODITB: Table where the terms of constraints are added */
2737 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2738 /* with ui and vj positive roots of the polynom of Legendre */
2739 /* of degree NBPNTU and NBPNTV respectively. */
2740 /* DIDITB: Table where the terms of constraints are added */
2741 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2742 /* with ui and vj positive roots of the polynom of Legendre */
2743 /* of degree NBPNTU and NBPNTV respectively. */
2745 /* COMMONS USED : */
2746 /* ---------------- */
2748 /* REFERENCES CALLED : */
2749 /* ----------------------- */
2751 /* DESCRIPTION/NOTES/LIMITATIONS : */
2752 /* ----------------------------------- */
2756 /* **********************************************************************
2759 /* Name of the routine */
2762 /* Parameter adjustments */
2763 diditb_dim1 = *nbpntu / 2 + 1;
2764 diditb_dim2 = *nbpntv / 2 + 1;
2765 diditb_offset = diditb_dim1 * diditb_dim2;
2766 diditb -= diditb_offset;
2767 disotb_dim1 = *nbpntu / 2;
2768 disotb_dim2 = *nbpntv / 2;
2769 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2770 disotb -= disotb_offset;
2771 soditb_dim1 = *nbpntu / 2;
2772 soditb_dim2 = *nbpntv / 2;
2773 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2774 soditb -= soditb_offset;
2775 sosotb_dim1 = *nbpntu / 2 + 1;
2776 sosotb_dim2 = *nbpntv / 2 + 1;
2777 sosotb_offset = sosotb_dim1 * sosotb_dim2;
2778 sosotb -= sosotb_offset;
2780 vhermt_dim1 = (*iordrv << 1) + 2;
2781 vhermt_offset = vhermt_dim1;
2782 vhermt -= vhermt_offset;
2783 fpntab_dim1 = *nbpntv;
2784 fpntab_offset = fpntab_dim1 + 1;
2785 fpntab -= fpntab_offset;
2786 ditbv2_dim1 = *nbpntu / 2 + 1;
2787 ditbv2_dim2 = *ndimen;
2788 ditbv2_offset = ditbv2_dim1 * (ditbv2_dim2 + 1);
2789 ditbv2 -= ditbv2_offset;
2790 ditbv1_dim1 = *nbpntu / 2 + 1;
2791 ditbv1_dim2 = *ndimen;
2792 ditbv1_offset = ditbv1_dim1 * (ditbv1_dim2 + 1);
2793 ditbv1 -= ditbv1_offset;
2794 sotbv2_dim1 = *nbpntu / 2 + 1;
2795 sotbv2_dim2 = *ndimen;
2796 sotbv2_offset = sotbv2_dim1 * (sotbv2_dim2 + 1);
2797 sotbv2 -= sotbv2_offset;
2798 sotbv1_dim1 = *nbpntu / 2 + 1;
2799 sotbv1_dim2 = *ndimen;
2800 sotbv1_offset = sotbv1_dim1 * (sotbv1_dim2 + 1);
2801 sotbv1 -= sotbv1_offset;
2804 ibb = AdvApp2Var_SysBase::mnfndeb_();
2806 AdvApp2Var_SysBase::mgenmsg_("MMA2CD2", 7L);
2809 /* ------------------- Discretization of Hermit polynoms -----------
2812 ncfhv = (*iordrv + 1) << 1;
2814 for (ii = 1; ii <= i__1; ++ii) {
2816 for (jj = 1; jj <= i__2; ++jj) {
2817 AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[ii * vhermt_dim1], &
2818 vrootl[jj], &fpntab[jj + ii * fpntab_dim1]);
2824 /* ---- The discretizations of polynoms of constraints are subtracted ----
2827 nuroo = *nbpntu / 2;
2828 nvroo = *nbpntv / 2;
2831 for (nd = 1; nd <= i__1; ++nd) {
2833 for (ii = 1; ii <= i__2; ++ii) {
2836 for (kk = 1; kk <= i__3; ++kk) {
2837 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1];
2838 bid2 = sotbv2[kk + (nd + ii * sotbv2_dim2) * sotbv2_dim1];
2839 bid3 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1];
2840 bid4 = ditbv2[kk + (nd + ii * ditbv2_dim2) * ditbv2_dim1];
2842 for (jj = 1; jj <= i__4; ++jj) {
2843 jjp = (*nbpntv + 1) / 2 + jj;
2844 jjm = nvroo - jj + 1;
2845 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] =
2846 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
2847 - bid1 * (fpntab[jjp + ((ii << 1) - 1) *
2848 fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) *
2849 fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) *
2850 fpntab_dim1] + fpntab[jjm + (ii << 1) *
2852 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] =
2853 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
2854 - bid3 * (fpntab[jjp + ((ii << 1) - 1) *
2855 fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) *
2856 fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) *
2857 fpntab_dim1] + fpntab[jjm + (ii << 1) *
2859 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] =
2860 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
2861 - bid1 * (fpntab[jjp + ((ii << 1) - 1) *
2862 fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) *
2863 fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) *
2864 fpntab_dim1] - fpntab[jjm + (ii << 1) *
2866 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] =
2867 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
2868 - bid3 * (fpntab[jjp + ((ii << 1) - 1) *
2869 fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) *
2870 fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) *
2871 fpntab_dim1] - fpntab[jjm + (ii << 1) *
2880 /* ------------ Case when the discretization is done only on the roots */
2881 /* ---------- of Legendre polynom of uneven degree, 0 is root */
2884 if (*nbpntv % 2 == 1) {
2886 for (ii = 1; ii <= i__2; ++ii) {
2888 for (kk = 1; kk <= i__3; ++kk) {
2889 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1]
2890 * fpntab[nvroo + 1 + ((ii << 1) - 1) *
2891 fpntab_dim1] + sotbv2[kk + (nd + ii * sotbv2_dim2)
2892 * sotbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) *
2894 sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2895 bid2 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1]
2896 * fpntab[nvroo + 1 + ((ii << 1) - 1) *
2897 fpntab_dim1] + ditbv2[kk + (nd + ii * ditbv2_dim2)
2898 * ditbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) *
2900 diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
2907 if (*nbpntu % 2 == 1) {
2909 for (ii = 1; ii <= i__2; ++ii) {
2911 for (jj = 1; jj <= i__3; ++jj) {
2912 jjp = (*nbpntv + 1) / 2 + jj;
2913 jjm = nvroo - jj + 1;
2914 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2915 fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] +
2916 fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) +
2917 sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2918 fpntab[jjp + (ii << 1) * fpntab_dim1] + fpntab[
2919 jjm + (ii << 1) * fpntab_dim1]);
2920 sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
2921 bid2 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2922 fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] -
2923 fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) +
2924 sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2925 fpntab[jjp + (ii << 1) * fpntab_dim1] - fpntab[
2926 jjm + (ii << 1) * fpntab_dim1]);
2927 diditb[jj + nd * diditb_dim2 * diditb_dim1] -= bid2;
2934 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2936 for (ii = 1; ii <= i__2; ++ii) {
2937 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * fpntab[
2938 nvroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbv2[(
2939 nd + ii * sotbv2_dim2) * sotbv2_dim1] * fpntab[nvroo
2940 + 1 + (ii << 1) * fpntab_dim1];
2941 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2950 /* ------------------------------ The End -------------------------------
2955 AdvApp2Var_SysBase::mgsomsg_("MMA2CD2", 7L);
2960 //=======================================================================
2961 //function : mma2cd3_
2963 //=======================================================================
2964 int mma2cd3_(integer *ndimen,
2983 /* System generated locals */
2984 integer sotbu1_dim1, sotbu1_dim2, sotbu1_offset, sotbu2_dim1, sotbu2_dim2,
2985 sotbu2_offset, ditbu1_dim1, ditbu1_dim2, ditbu1_offset,
2986 ditbu2_dim1, ditbu2_dim2, ditbu2_offset, fpntab_dim1,
2987 fpntab_offset, uhermt_dim1, uhermt_offset, sosotb_dim1,
2988 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2989 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2990 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2992 /* Local variables */
2993 integer ncfhu, nuroo, nvroo, ii, nd, jj, kk, ibb, kkm, kkp;
2994 doublereal bid1, bid2, bid3, bid4;
2996 /* **********************************************************************
3000 /* Discretisation on the parameters of polynoms of interpolation */
3001 /* of constraints on 2 borders iso-U of order IORDRU. */
3006 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3008 /* INPUT ARGUMENTS : */
3009 /* ------------------ */
3010 /* NDIMEN: Dimension of the space. */
3011 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3012 /* This is also the nb of root of Legendre polynom where discretization is done. */
3013 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
3015 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
3016 /* This is also the nb of root of Legendre polynom where discretization is done. */
3017 /* IORDRV: Order of constraint imposed at the extremities of iso-V */
3018 /* = 0, calculate the extremities of iso-V */
3019 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3020 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3021 /* SOTBU1: Table of NBPNTU/2 sums of 2 index points */
3022 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3023 /* SOTBU2: Table of NBPNTV/2 sums of 2 index points */
3024 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3025 /* DITBU1: Table of NBPNTU/2 differences of 2 index points */
3026 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3027 /* DITBU2: Table of NBPNTU/2 differences of 2 index points */
3028 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3029 /* SOSOTB: Preinitialized table (input/output argument). */
3030 /* DISOTB: Preinitialized table (input/output argument). */
3031 /* SODITB: Preinitialized table (input/output argument). */
3032 /* DIDITB: Preinitialized table (input/output argument) */
3034 /* OUTPUT ARGUMENTS : */
3035 /* ------------------- */
3036 /* FPNTAB: Auxiliary table. */
3037 /* UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
3038 /* SOSOTB: Table where the terms of constraints are added */
3039 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3040 /* with ui and vj positive roots of the Legendre polynom */
3041 /* of degree NBPNTU and NBPNTV respectively. */
3042 /* DISOTB: Table where the terms of constraints are added */
3043 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3044 /* with ui and vj positive roots of the polynom of Legendre */
3045 /* of degree NBPNTU and NBPNTV respectively. */
3046 /* SODITB: Table where the terms of constraints are added */
3047 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3048 /* with ui and vj positive roots of the polynom of Legendre */
3049 /* of degree NBPNTU and NBPNTV respectively. */
3050 /* DIDITB: Table where the terms of constraints are added */
3051 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3052 /* with ui and vj positive roots of the polynom of Legendre */
3053 /* of degree NBPNTU and NBPNTV respectively. */
3055 /* COMMONS USED : */
3056 /* ---------------- */
3058 /* REFERENCES CALLED : */
3059 /* ----------------------- */
3061 /* DESCRIPTION/NOTES/LIMITATIONS : */
3062 /* ----------------------------------- */
3064 /* $ HISTORIQUE DES MODIFICATIONS : */
3065 /* -------------------------------- */
3066 /* 08-08-1991: RBD; Creation. */
3068 /* **********************************************************************
3071 /* Name of the routine */
3074 /* Parameter adjustments */
3076 diditb_dim1 = *nbpntu / 2 + 1;
3077 diditb_dim2 = *nbpntv / 2 + 1;
3078 diditb_offset = diditb_dim1 * diditb_dim2;
3079 diditb -= diditb_offset;
3080 disotb_dim1 = *nbpntu / 2;
3081 disotb_dim2 = *nbpntv / 2;
3082 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3083 disotb -= disotb_offset;
3084 soditb_dim1 = *nbpntu / 2;
3085 soditb_dim2 = *nbpntv / 2;
3086 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3087 soditb -= soditb_offset;
3088 sosotb_dim1 = *nbpntu / 2 + 1;
3089 sosotb_dim2 = *nbpntv / 2 + 1;
3090 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3091 sosotb -= sosotb_offset;
3092 uhermt_dim1 = (*iordru << 1) + 2;
3093 uhermt_offset = uhermt_dim1;
3094 uhermt -= uhermt_offset;
3095 fpntab_dim1 = *nbpntu;
3096 fpntab_offset = fpntab_dim1 + 1;
3097 fpntab -= fpntab_offset;
3098 ditbu2_dim1 = *nbpntv / 2 + 1;
3099 ditbu2_dim2 = *ndimen;
3100 ditbu2_offset = ditbu2_dim1 * (ditbu2_dim2 + 1);
3101 ditbu2 -= ditbu2_offset;
3102 ditbu1_dim1 = *nbpntv / 2 + 1;
3103 ditbu1_dim2 = *ndimen;
3104 ditbu1_offset = ditbu1_dim1 * (ditbu1_dim2 + 1);
3105 ditbu1 -= ditbu1_offset;
3106 sotbu2_dim1 = *nbpntv / 2 + 1;
3107 sotbu2_dim2 = *ndimen;
3108 sotbu2_offset = sotbu2_dim1 * (sotbu2_dim2 + 1);
3109 sotbu2 -= sotbu2_offset;
3110 sotbu1_dim1 = *nbpntv / 2 + 1;
3111 sotbu1_dim2 = *ndimen;
3112 sotbu1_offset = sotbu1_dim1 * (sotbu1_dim2 + 1);
3113 sotbu1 -= sotbu1_offset;
3116 ibb = AdvApp2Var_SysBase::mnfndeb_();
3118 AdvApp2Var_SysBase::mgenmsg_("MMA2CD3", 7L);
3121 /* ------------------- Discretization of polynoms of Hermit -----------
3124 ncfhu = (*iordru + 1) << 1;
3126 for (ii = 1; ii <= i__1; ++ii) {
3128 for (kk = 1; kk <= i__2; ++kk) {
3129 AdvApp2Var_MathBase::mmmpocur_(&ncfhu,
3132 &uhermt[ii * uhermt_dim1],
3134 &fpntab[kk + ii * fpntab_dim1]);
3140 /* ---- The discretizations of polynoms of constraints are subtracted ----
3143 nvroo = *nbpntv / 2;
3144 nuroo = *nbpntu / 2;
3147 for (nd = 1; nd <= i__1; ++nd) {
3149 for (ii = 1; ii <= i__2; ++ii) {
3152 for (jj = 1; jj <= i__3; ++jj) {
3153 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1];
3154 bid2 = sotbu2[jj + (nd + ii * sotbu2_dim2) * sotbu2_dim1];
3155 bid3 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1];
3156 bid4 = ditbu2[jj + (nd + ii * ditbu2_dim2) * ditbu2_dim1];
3158 for (kk = 1; kk <= i__4; ++kk) {
3159 kkp = (*nbpntu + 1) / 2 + kk;
3160 kkm = nuroo - kk + 1;
3161 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] =
3162 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
3163 - bid1 * (fpntab[kkp + ((ii << 1) - 1) *
3164 fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) *
3165 fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) *
3166 fpntab_dim1] + fpntab[kkm + (ii << 1) *
3168 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] =
3169 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
3170 - bid1 * (fpntab[kkp + ((ii << 1) - 1) *
3171 fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) *
3172 fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) *
3173 fpntab_dim1] - fpntab[kkm + (ii << 1) *
3175 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] =
3176 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
3177 - bid3 * (fpntab[kkp + ((ii << 1) - 1) *
3178 fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) *
3179 fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) *
3180 fpntab_dim1] + fpntab[kkm + (ii << 1) *
3182 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] =
3183 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
3184 - bid3 * (fpntab[kkp + ((ii << 1) - 1) *
3185 fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) *
3186 fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) *
3187 fpntab_dim1] - fpntab[kkm + (ii << 1) *
3196 /* ------------ Case when the discretization is done only on the roots */
3197 /* ---------- of Legendre polynom of uneven degree, 0 is root */
3201 if (*nbpntu % 2 == 1) {
3203 for (ii = 1; ii <= i__2; ++ii) {
3205 for (jj = 1; jj <= i__3; ++jj) {
3206 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1]
3207 * fpntab[nuroo + 1 + ((ii << 1) - 1) *
3208 fpntab_dim1] + sotbu2[jj + (nd + ii * sotbu2_dim2)
3209 * sotbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) *
3211 sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
3212 bid2 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1]
3213 * fpntab[nuroo + 1 + ((ii << 1) - 1) *
3214 fpntab_dim1] + ditbu2[jj + (nd + ii * ditbu2_dim2)
3215 * ditbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) *
3217 diditb[(jj + nd * diditb_dim2) * diditb_dim1] -= bid2;
3224 if (*nbpntv % 2 == 1) {
3226 for (ii = 1; ii <= i__2; ++ii) {
3228 for (kk = 1; kk <= i__3; ++kk) {
3229 kkp = (*nbpntu + 1) / 2 + kk;
3230 kkm = nuroo - kk + 1;
3231 bid1 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * (
3232 fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] +
3233 fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1]) +
3234 sotbu2[(nd + ii * sotbu2_dim2) * sotbu2_dim1] * (
3235 fpntab[kkp + (ii << 1) * fpntab_dim1] + fpntab[
3236 kkm + (ii << 1) * fpntab_dim1]);
3237 sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
3238 bid2 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * (
3239 fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] -
3240 fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1]) +
3241 sotbu2[(nd + ii * sotbu2_dim2) * sotbu2_dim1] * (
3242 fpntab[kkp + (ii << 1) * fpntab_dim1] - fpntab[
3243 kkm + (ii << 1) * fpntab_dim1]);
3244 diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
3251 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
3253 for (ii = 1; ii <= i__2; ++ii) {
3254 bid1 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * fpntab[
3255 nuroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbu2[(
3256 nd + ii * sotbu2_dim2) * sotbu2_dim1] * fpntab[nuroo
3257 + 1 + (ii << 1) * fpntab_dim1];
3258 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
3267 /* ------------------------------ The End -------------------------------
3272 AdvApp2Var_SysBase::mgsomsg_("MMA2CD3", 7L);
3277 //=======================================================================
3278 //function : mma2cdi_
3280 //=======================================================================
3281 int AdvApp2Var_ApproxF2var::mma2cdi_( integer *ndimen,
3309 /* System generated locals */
3310 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
3311 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
3312 contr4_dim1, contr4_dim2, contr4_offset, sosotb_dim1, sosotb_dim2,
3313 sosotb_offset, diditb_dim1, diditb_dim2, diditb_offset,
3314 soditb_dim1, soditb_dim2, soditb_offset, disotb_dim1, disotb_dim2,
3317 /* Local variables */
3320 doublereal* wrkar = 0;
3322 integer ibb, ier = 0;
3323 integer isz1, isz2, isz3, isz4;
3324 intptr_t ipt1, ipt2, ipt3, ipt4;
3329 /* **********************************************************************
3334 /* Discretisation on the parameters of polynomes of interpolation */
3335 /* of constraints of order IORDRE. */
3339 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3341 //* INPUT ARGUMENTS : */
3342 /* ------------------ */
3343 /* NDIMEN: Dimension of the space. */
3344 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3345 /* This is also the nb of root of Legendre polynom where discretization is done. */
3346 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
3348 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
3349 /* This is also the nb of root of Legendre polynom where discretization is done. */
3350 /* VROOTL: Table of parameters of discretisation ON (-1,1) by V.*/
3352 /* IORDRV: Order of constraint imposed at the extremities of iso-U */
3353 /* = 0, calculate the extremities of iso-U */
3354 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
3355 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
3356 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
3357 /* = 0, calculate the extremities of iso-V */
3358 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3359 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3360 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
3361 /* extremities of F(U0,V0) and its derivatives. */
3362 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
3363 /* extremities of F(U1,V0) and its derivatives. */
3364 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
3365 /* extremities of F(U0,V1) and its derivatives. */
3366 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
3367 /* extremities of F(U1,V1) and its derivatives. */
3368 /* SOTBU1: Table of NBPNTU/2 sums of 2 index points */
3369 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3370 /* SOTBU2: Table of NBPNTV/2 sums of 2 index points */
3371 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3372 /* DITBU1: Table of NBPNTU/2 differences of 2 index points */
3373 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3374 /* DITBU2: Table of NBPNTU/2 differences of 2 index points */
3375 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3376 /* SOTBV1: Table of NBPNTV/2 sums of 2 index points */
3377 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
3378 /* SOTBV2: Table of NBPNTV/2 sums of 2 index points */
3379 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
3380 /* DITBV1: Table of NBPNTV/2 differences of 2 index points */
3381 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
3382 /* DITBV2: Table of NBPNTV/2 differences of 2 index points */
3383 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
3384 /* SOSOTB: Preinitialized table (input/output argument). */
3385 /* DISOTB: Preinitialized table (input/output argument). */
3386 /* SODITB: Preinitialized table (input/output argument). */
3387 /* DIDITB: Preinitialized table (input/output argument) */
3389 /* ARGUMENTS DE SORTIE : */
3390 /* ------------------- */
3391 /* SOSOTB: Table where the terms of constraints are added */
3392 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3393 /* with ui and vj positive roots of the Legendre polynom */
3394 /* of degree NBPNTU and NBPNTV respectively. */
3395 /* DISOTB: Table where the terms of constraints are added */
3396 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3397 /* with ui and vj positive roots of the polynom of Legendre */
3398 /* of degree NBPNTU and NBPNTV respectively. */
3399 /* SODITB: Table where the terms of constraints are added */
3400 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3401 /* with ui and vj positive roots of the polynom of Legendre */
3402 /* of degree NBPNTU and NBPNTV respectively. */
3403 /* DIDITB: Table where the terms of constraints are added */
3404 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3405 /* with ui and vj positive roots of the polynom of Legendre */
3406 /* of degree NBPNTU and NBPNTV respectively. */
3407 /* IERCOD: = 0, OK, */
3408 /* = 1, Value or IORDRV or IORDRU is out of allowed values. */
3409 /* =13, Pb of dynamic allocation. */
3411 /* COMMONS USED : */
3412 /* ---------------- */
3414 /* REFERENCES CALLED : */
3415 /* -------------------- */
3417 /* DESCRIPTION/NOTES/LIMITATIONS : */
3418 /* ------------------------------- */
3421 /* **********************************************************************
3424 /* The name of the routine */
3427 /* Parameter adjustments */
3429 diditb_dim1 = *nbpntu / 2 + 1;
3430 diditb_dim2 = *nbpntv / 2 + 1;
3431 diditb_offset = diditb_dim1 * diditb_dim2;
3432 diditb -= diditb_offset;
3433 disotb_dim1 = *nbpntu / 2;
3434 disotb_dim2 = *nbpntv / 2;
3435 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3436 disotb -= disotb_offset;
3437 soditb_dim1 = *nbpntu / 2;
3438 soditb_dim2 = *nbpntv / 2;
3439 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3440 soditb -= soditb_offset;
3441 sosotb_dim1 = *nbpntu / 2 + 1;
3442 sosotb_dim2 = *nbpntv / 2 + 1;
3443 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3444 sosotb -= sosotb_offset;
3446 contr4_dim1 = *ndimen;
3447 contr4_dim2 = *iordru + 2;
3448 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
3449 contr4 -= contr4_offset;
3450 contr3_dim1 = *ndimen;
3451 contr3_dim2 = *iordru + 2;
3452 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
3453 contr3 -= contr3_offset;
3454 contr2_dim1 = *ndimen;
3455 contr2_dim2 = *iordru + 2;
3456 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
3457 contr2 -= contr2_offset;
3458 contr1_dim1 = *ndimen;
3459 contr1_dim2 = *iordru + 2;
3460 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
3461 contr1 -= contr1_offset;
3470 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
3473 ibb = AdvApp2Var_SysBase::mnfndeb_();
3475 AdvApp2Var_SysBase::mgenmsg_("MMA2CDI", 7L);
3479 if (*iordru < -1 || *iordru > 2) {
3482 if (*iordrv < -1 || *iordrv > 2) {
3486 /* ------------------------- Set to zero --------------------------------
3489 ilong = (*nbpntu / 2 + 1) * (*nbpntv / 2 + 1) * *ndimen;
3490 AdvApp2Var_SysBase::mvriraz_(&ilong, &sosotb[sosotb_offset]);
3491 AdvApp2Var_SysBase::mvriraz_(&ilong, &diditb[diditb_offset]);
3492 ilong = *nbpntu / 2 * (*nbpntv / 2) * *ndimen;
3493 AdvApp2Var_SysBase::mvriraz_(&ilong, &soditb[soditb_offset]);
3494 AdvApp2Var_SysBase::mvriraz_(&ilong, &disotb[disotb_offset]);
3495 if (*iordru == -1 && *iordrv == -1) {
3501 isz1 = ((*iordru + 1) << 2) * (*iordru + 1);
3502 isz2 = ((*iordrv + 1) << 2) * (*iordrv + 1);
3503 isz3 = ((*iordru + 1) << 1) * *nbpntu;
3504 isz4 = ((*iordrv + 1) << 1) * *nbpntv;
3505 iszwr = isz1 + isz2 + isz3 + isz4;
3506 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3515 if (*iordru >= 0 && *iordru <= 2) {
3517 /* --- Return 2*(IORDRU+1) coeff of 2*(IORDRU+1) polynoms of Hermite
3520 AdvApp2Var_ApproxF2var::mma1her_(iordru, &wrkar[ipt1], iercod);
3525 /* ---- Subract discretizations of polynoms of constraints
3528 mma2cd3_(ndimen, nbpntu, &urootl[1], nbpntv, iordru, &sotbu1[1], &
3529 sotbu2[1], &ditbu1[1], &ditbu2[1], &wrkar[ipt3], &wrkar[ipt1],
3530 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3531 disotb_offset], &diditb[diditb_offset]);
3534 if (*iordrv >= 0 && *iordrv <= 2) {
3536 /* --- Return 2*(IORDRV+1) coeff of 2*(IORDRV+1) polynoms of Hermite
3539 AdvApp2Var_ApproxF2var::mma1her_(iordrv, &wrkar[ipt2], iercod);
3544 /* ---- Subtract discretisations of polynoms of constraint
3547 mma2cd2_(ndimen, nbpntu, nbpntv, &vrootl[1], iordrv, &sotbv1[1], &
3548 sotbv2[1], &ditbv1[1], &ditbv2[1], &wrkar[ipt4], &wrkar[ipt2],
3549 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3550 disotb_offset], &diditb[diditb_offset]);
3553 /* --------------- Subtract constraints of corners ----------------
3556 if (*iordru >= 0 && *iordrv >= 0) {
3557 mma2cd1_(ndimen, nbpntu, &urootl[1], nbpntv, &vrootl[1], iordru,
3558 iordrv, &contr1[contr1_offset], &contr2[contr2_offset], &
3559 contr3[contr3_offset], &contr4[contr4_offset], &wrkar[ipt3], &
3560 wrkar[ipt4], &wrkar[ipt1], &wrkar[ipt2], &sosotb[
3561 sosotb_offset], &soditb[soditb_offset], &disotb[disotb_offset]
3562 , &diditb[diditb_offset]);
3566 /* ------------------------------ The End -------------------------------
3568 /* --> IORDRE is not within the autorised diapason. */
3572 /* --> PB of dynamic allocation. */
3579 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3584 AdvApp2Var_SysBase::maermsg_("MMA2CDI", iercod, 7L);
3586 AdvApp2Var_SysBase::mgsomsg_("MMA2CDI", 7L);
3591 //=======================================================================
3592 //function : mma2ce1_
3594 //=======================================================================
3595 int AdvApp2Var_ApproxF2var::mma2ce1_(integer *numdec,
3625 /* System generated locals */
3626 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3627 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3628 diditb_dim1, diditb_dim2, diditb_offset, patjac_dim1, patjac_dim2,
3631 /* Local variables */
3634 doublereal* wrkar = 0;
3637 integer isz1, isz2, isz3, isz4, isz5, isz6, isz7;
3638 intptr_t ipt1, ipt2, ipt3, ipt4, ipt5, ipt6, ipt7;
3642 /* **********************************************************************
3647 /* Calculation of coefficients of polynomial approximation of degree */
3648 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3649 /* discretization on roots of Legendre polynom of degree */
3650 /* NBPNTU by U and NBPNTV by V. */
3654 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&POLYNOME,&ERREUR */
3656 /* INPUT ARGUMENTS : */
3657 /* ------------------ */
3658 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3659 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3660 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3661 /* directions simultaneously (cutting by V is preferable). */
3662 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3663 /* directions simultaneously (cutting by U is preferable). */
3664 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3665 /* of cutting Vj). */
3666 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3667 /* of cutting Ui). */
3668 /* = 0, It is not POSSIBLE to cut anything */
3669 /* NDIMEN: Dimension of the space. */
3670 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3671 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3672 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3673 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3674 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3675 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3676 /* NDJACU: Max degree of the polynom of approximation by U. */
3677 /* The representation in the orthogonal base starts from degree */
3678 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3679 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3680 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3681 /* NDJACV: Max degree of the polynom of approximation by V. */
3682 /* The representation in the orthogonal base starts from degree */
3683 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3684 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3685 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3686 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3687 /* to the step of constraints C0, C1 or C2. */
3688 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3689 /* to the step of constraints C0, C1 or C2. */
3690 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3691 /* calculated the coefficients of integration by u */
3692 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3693 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3694 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3695 /* calculated the coefficients of integration by u */
3696 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3697 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3698 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3699 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3700 /* with ui and vj - positive roots of the Legendre polynom */
3701 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3702 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3703 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3704 /* SOSOTB(0,0) contains F(0,0). */
3705 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3706 /* with ui and vj positive roots of Legendre polynom */
3707 /* of degree NBPNTU and NBPNTV respectively. */
3708 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3709 /* with ui and vj positive roots of Legendre polynom */
3710 /* of degree NBPNTU and NBPNTV respectively. */
3711 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3712 /* with ui and vj positive roots of Legendre polynom */
3713 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3714 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3715 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3717 /* OUTPUT ARGUMENTS */
3718 /* --------------- */
3719 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
3720 /* of F(u,v) with eventually taking into account of */
3721 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
3722 /* This table contains other coeff if ITYDEC = 0. */
3723 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
3724 /* on each of sub-spaces SI ITYDEC = 0. */
3725 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
3726 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
3727 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
3728 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
3729 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
3730 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
3731 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
3732 /* = 3, it is NECESSARY to cut both by U AND by V. */
3733 /* IERCOD: Error code. */
3734 /* = 0, Everything is OK. */
3735 /* = -1, There is the best possible solution, but the */
3736 /* user tolerance is not satisfactory (3*only) */
3737 /* = 1, Incoherent entries. */
3739 /* COMMONS USED : */
3740 /* ---------------- */
3742 /* REFERENCES CALLED : */
3743 /* --------------------- */
3745 /* DESCRIPTION/NOTES/LIMITATIONS : */
3746 /* ------------------------------- */
3749 /* **********************************************************************
3751 /* Name of the routine */
3754 /* --------------------------- Initialisations --------------------------
3757 /* Parameter adjustments */
3762 patjac_dim1 = *ndjacu + 1;
3763 patjac_dim2 = *ndjacv + 1;
3764 patjac_offset = patjac_dim1 * patjac_dim2;
3765 patjac -= patjac_offset;
3766 diditb_dim1 = *nbpntu / 2 + 1;
3767 diditb_dim2 = *nbpntv / 2 + 1;
3768 diditb_offset = diditb_dim1 * diditb_dim2;
3769 diditb -= diditb_offset;
3770 soditb_dim1 = *nbpntu / 2;
3771 soditb_dim2 = *nbpntv / 2;
3772 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3773 soditb -= soditb_offset;
3774 disotb_dim1 = *nbpntu / 2;
3775 disotb_dim2 = *nbpntv / 2;
3776 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3777 disotb -= disotb_offset;
3778 sosotb_dim1 = *nbpntu / 2 + 1;
3779 sosotb_dim2 = *nbpntv / 2 + 1;
3780 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3781 sosotb -= sosotb_offset;
3784 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
3786 AdvApp2Var_SysBase::mgenmsg_("MMA2CE1", 7L);
3791 isz1 = (*nbpntu / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1);
3792 isz2 = (*nbpntv / 2 + 1) * (*ndjacv - ((*iordrv + 1) << 1) + 1);
3793 isz3 = (*nbpntv / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3794 isz4 = *nbpntv / 2 * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3795 isz5 = *ndjacu + 1 - ((*iordru + 1) << 1);
3796 isz6 = *ndjacv + 1 - ((*iordrv + 1) << 1);
3797 isz7 = *ndimen << 2;
3798 iszwr = isz1 + isz2 + isz3 + isz4 + isz5 + isz6 + isz7;
3799 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
3800 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3812 /* ----------------- Return Gauss coefficients of integration ----------------
3815 AdvApp2Var_ApproxF2var::mmapptt_(ndjacu, nbpntu, iordru, &wrkar[ipt1], iercod);
3819 AdvApp2Var_ApproxF2var::mmapptt_(ndjacv, nbpntv, iordrv, &wrkar[ipt2], iercod);
3824 /* ------------------- Return max polynoms of Jacobi ------------
3827 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacu, iordru, &wrkar[ipt5]);
3828 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacv, iordrv, &wrkar[ipt6]);
3830 /* ------ Calculate the coefficients and their contribution to the error ----
3833 mma2ce2_(numdec, ndimen, nbsesp, &ndimse[1], ndminu, ndminv, ndguli,
3834 ndgvli, ndjacu, ndjacv, iordru, iordrv, nbpntu, nbpntv, &epsapr[1]
3835 , &sosotb[sosotb_offset], &disotb[disotb_offset], &soditb[
3836 soditb_offset], &diditb[diditb_offset], &wrkar[ipt1], &wrkar[ipt2]
3837 , &wrkar[ipt5], &wrkar[ipt6], &wrkar[ipt7], &wrkar[ipt3], &wrkar[
3838 ipt4], &patjac[patjac_offset], &errmax[1], &errmoy[1], ndegpu,
3839 ndegpv, itydec, iercod);
3845 /* ------------------------------ The end -------------------------------
3854 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3859 AdvApp2Var_SysBase::maermsg_("MMA2CE1", iercod, 7L);
3861 AdvApp2Var_SysBase::mgsomsg_("MMA2CE1", 7L);
3866 //=======================================================================
3867 //function : mma2ce2_
3869 //=======================================================================
3870 int mma2ce2_(integer *numdec,
3905 /* System generated locals */
3906 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3907 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3908 diditb_dim1, diditb_dim2, diditb_offset, gssutb_dim1, gssvtb_dim1,
3909 chpair_dim1, chpair_dim2, chpair_offset, chimpr_dim1,
3910 chimpr_dim2, chimpr_offset, patjac_dim1, patjac_dim2,
3911 patjac_offset, vecerr_dim1, vecerr_offset, i__1, i__2, i__3, i__4;
3913 /* Local variables */
3915 integer idim, igsu, minu, minv, maxu, maxv, igsv;
3917 integer i2rdu, i2rdv, ndses, nd, ii, jj, kk, nu, nv;
3921 /* **********************************************************************
3925 /* Calculation of coefficients of polynomial approximation of degree */
3926 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3927 /* discretization on roots of Legendre polynom of degree */
3928 /* NBPNTU by U and NBPNTV by V. */
3932 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&COEFFICIENT,&POLYNOME */
3934 /* INPUT ARGUMENTS : */
3935 /* ------------------ */
3936 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3937 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3938 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3939 /* directions simultaneously (cutting by V is preferable). */
3940 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3941 /* directions simultaneously (cutting by U is preferable). */
3942 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3943 /* of cutting Vj). */
3944 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3945 /* of cutting Ui). */
3946 /* = 0, It is not POSSIBLE to cut anything */
3947 /* NDIMEN: Total dimension of the space. */
3948 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3949 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3950 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3951 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3952 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3953 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3954 /* NDJACU: Max degree of the polynom of approximation by U. */
3955 /* The representation in the orthogonal base starts from degree */
3956 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3957 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3958 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3959 /* NDJACV: Max degree of the polynom of approximation by V. */
3960 /* The representation in the orthogonal base starts from degree */
3961 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3962 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3963 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3964 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3965 /* to the step of constraints C0, C1 or C2. */
3966 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3967 /* to the step of constraints C0, C1 or C2. */
3968 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3969 /* calculated the coefficients of integration by u */
3970 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3971 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3972 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3973 /* calculated the coefficients of integration by u */
3974 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3975 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3976 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3977 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3978 /* with ui and vj - positive roots of the Legendre polynom */
3979 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3980 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3981 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3982 /* SOSOTB(0,0) contains F(0,0). */
3983 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3984 /* with ui and vj positive roots of Legendre polynom */
3985 /* of degree NBPNTU and NBPNTV respectively. */
3986 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3987 /* with ui and vj positive roots of Legendre polynom */
3988 /* of degree NBPNTU and NBPNTV respectively. */
3989 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3990 /* with ui and vj positive roots of Legendre polynom */
3991 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3992 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3993 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3994 /* GSSUTB: Table of coefficients of integration by Gauss method */
3995 /* by U: i varies from 0 to NBPNTU/2 and k varies from 0 to */
3996 /* NDJACU-2*(IORDRU+1). */
3997 /* GSSVTB: Table of coefficients of integration by Gauss method */
3998 /* by V: i varies from 0 to NBPNTV/2 and k varies from 0 to */
3999 /* NDJACV-2*(IORDRV+1). */
4000 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
4001 /* from degree 0 to degree NDJACU - 2*(IORDRU+1) */
4002 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
4003 /* from degree 0 to degree NDJACV - 2*(IORDRV+1) */
4005 /* OUTPUT ARGUMENTS : */
4006 /* ------------------- */
4007 /* VECERR: Auxiliary table. */
4008 /* CHPAIR: Auxiliary table of terms connected to degree NDJACU by U */
4009 /* to calculate the coeff. of approximation of EVEN degree by V. */
4010 /* CHIMPR: Auxiliary table of terms connected to degree NDJACU by U */
4011 /* to calculate the coeff. of approximation of UNEVEN degree by V. */
4012 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
4013 /* of F(u,v) with eventually taking into account of */
4014 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
4015 /* This table contains other coeff if ITYDEC = 0. */
4016 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
4017 /* on each of sub-spaces SI ITYDEC = 0. */
4018 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
4019 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
4020 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
4021 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
4022 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
4023 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
4024 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
4025 /* = 3, it is NECESSARY to cut both by U AND by V. */
4026 /* IERCOD: Error code. */
4027 /* = 0, Everything is OK. */
4028 /* = -1, There is the best possible solution, but the */
4029 /* user tolerance is not satisfactory (3*only) */
4030 /* = 1, Incoherent entries. */
4032 /* COMMONS USED : */
4033 /* ---------------- */
4035 /* REFERENCES CALLED : */
4036 /* --------------------- */
4038 /* DESCRIPTION/NOTES/LIMITATIONS : */
4040 /* **********************************************************************
4042 /* Name of the routine */
4045 /* --------------------------- Initialisations --------------------------
4048 /* Parameter adjustments */
4049 vecerr_dim1 = *ndimen;
4050 vecerr_offset = vecerr_dim1 + 1;
4051 vecerr -= vecerr_offset;
4056 patjac_dim1 = *ndjacu + 1;
4057 patjac_dim2 = *ndjacv + 1;
4058 patjac_offset = patjac_dim1 * patjac_dim2;
4059 patjac -= patjac_offset;
4060 gssutb_dim1 = *nbpntu / 2 + 1;
4061 chimpr_dim1 = *nbpntv / 2;
4062 chimpr_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4063 chimpr_offset = chimpr_dim1 * chimpr_dim2 + 1;
4064 chimpr -= chimpr_offset;
4065 chpair_dim1 = *nbpntv / 2 + 1;
4066 chpair_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4067 chpair_offset = chpair_dim1 * chpair_dim2;
4068 chpair -= chpair_offset;
4069 gssvtb_dim1 = *nbpntv / 2 + 1;
4070 diditb_dim1 = *nbpntu / 2 + 1;
4071 diditb_dim2 = *nbpntv / 2 + 1;
4072 diditb_offset = diditb_dim1 * diditb_dim2;
4073 diditb -= diditb_offset;
4074 soditb_dim1 = *nbpntu / 2;
4075 soditb_dim2 = *nbpntv / 2;
4076 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
4077 soditb -= soditb_offset;
4078 disotb_dim1 = *nbpntu / 2;
4079 disotb_dim2 = *nbpntv / 2;
4080 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
4081 disotb -= disotb_offset;
4082 sosotb_dim1 = *nbpntu / 2 + 1;
4083 sosotb_dim2 = *nbpntv / 2 + 1;
4084 sosotb_offset = sosotb_dim1 * sosotb_dim2;
4085 sosotb -= sosotb_offset;
4088 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4090 AdvApp2Var_SysBase::mgenmsg_("MMA2CE2", 7L);
4092 /* --> A priori everything is OK */
4094 /* --> test of inputs */
4095 if (*numdec < 0 || *numdec > 5) {
4098 if ((*iordru << 1) + 1 > *ndminu) {
4101 if (*ndminu > *ndguli) {
4104 if (*ndguli >= *ndjacu) {
4107 if ((*iordrv << 1) + 1 > *ndminv) {
4110 if (*ndminv > *ndgvli) {
4113 if (*ndgvli >= *ndjacv) {
4116 /* --> A priori, no cuts to be done */
4118 /* --> Min. degrees to return: NDMINU,NDMINV */
4121 /* --> For the moment, max errors are null */
4122 AdvApp2Var_SysBase::mvriraz_(nbsesp, &errmax[1]);
4124 AdvApp2Var_SysBase::mvriraz_(&nd, &vecerr[vecerr_offset]);
4125 /* --> and the square, too. */
4126 nd = (*ndjacu + 1) * (*ndjacv + 1) * *ndimen;
4127 AdvApp2Var_SysBase::mvriraz_(&nd, &patjac[patjac_offset]);
4129 i2rdu = (*iordru + 1) << 1;
4130 i2rdv = (*iordrv + 1) << 1;
4132 /* **********************************************************************
4134 /* -------------------- HERE IT IS POSSIBLE TO CUT ----------------------
4136 /* **********************************************************************
4139 if (*numdec > 0 && *numdec <= 5) {
4141 /* ******************************************************************
4143 /* ---------------------- Calculate coeff of zone 4 -------------
4157 /* ---------------- Calculate the terms connected to degree by U ---------
4161 for (nd = 1; nd <= i__1; ++nd) {
4163 for (kk = minu; kk <= i__2; ++kk) {
4165 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4166 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4167 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4168 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4169 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4170 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4171 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4177 /* ------------------- Calculate the coefficients of PATJAC ------------
4180 igsu = minu - i2rdu;
4182 for (jj = minv; jj <= i__1; ++jj) {
4185 for (nd = 1; nd <= i__2; ++nd) {
4186 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4187 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4188 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4189 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4190 patjac_dim2) * patjac_dim1]);
4194 /* ----- Contribution of calculated terms to the approximation error */
4195 /* for terms (I,J) with MINU <= I <= MAXU, J fixe. */
4199 for (nd = 1; nd <= i__2; ++nd) {
4201 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &jj, &jj,
4202 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4203 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4204 &vecerr[nd + (vecerr_dim1 << 2)]);
4205 if (vecerr[nd + (vecerr_dim1 << 2)] > epsapr[nd]) {
4214 /* ******************************************************************
4216 /* ---------------------- Calculate the coeff of zone 2 -------------
4219 minu = (*iordru + 1) << 1;
4224 /* --> If zone 2 is empty, pass to zone 3. */
4225 /* VECERR(ND,2) was already set to zero. */
4230 /* ---------------- Calculate the terms connected to degree by U ------------
4234 for (nd = 1; nd <= i__1; ++nd) {
4236 for (kk = minu; kk <= i__2; ++kk) {
4238 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4239 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4240 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4241 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4242 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4243 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4244 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4250 /* ------------------- Calculate the coefficients of PATJAC ------------
4253 igsu = minu - i2rdu;
4255 for (jj = minv; jj <= i__1; ++jj) {
4258 for (nd = 1; nd <= i__2; ++nd) {
4259 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4260 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4261 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4262 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4263 patjac_dim2) * patjac_dim1]);
4269 /* -----Contribution of calculated terms to the approximation error */
4270 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4274 for (nd = 1; nd <= i__1; ++nd) {
4276 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4277 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4278 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4279 vecerr[nd + (vecerr_dim1 << 1)]);
4284 /* ******************************************************************
4286 /* ---------------------- Calculation of coeff of zone 3 -------------
4292 minv = (*iordrv + 1) << 1;
4295 /* -> If zone 3 is empty, pass to the test of cutting. */
4296 /* VECERR(ND,3) was already set to zero */
4301 /* ----------- The terms connected to the degree by U are already calculated -----
4303 /* ------------------- Calculation of coefficients of PATJAC ------------
4306 igsu = minu - i2rdu;
4308 for (jj = minv; jj <= i__1; ++jj) {
4311 for (nd = 1; nd <= i__2; ++nd) {
4312 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4313 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4314 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4315 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4316 patjac_dim2) * patjac_dim1]);
4322 /* ----- Contribution of calculated terms to the approximation error */
4323 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV. */
4327 for (nd = 1; nd <= i__1; ++nd) {
4329 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4330 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4331 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4332 vecerr[nd + vecerr_dim1 * 3]);
4337 /* ******************************************************************
4339 /* --------------------------- Tests of cutting ---------------------
4344 for (nd = 1; nd <= i__1; ++nd) {
4345 vaux[0] = vecerr[nd + (vecerr_dim1 << 1)];
4346 vaux[1] = vecerr[nd + (vecerr_dim1 << 2)];
4347 vaux[2] = vecerr[nd + vecerr_dim1 * 3];
4349 errmax[nd] = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4350 if (errmax[nd] > epsapr[nd]) {
4352 zv = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4353 zu = AdvApp2Var_MathBase::mzsnorm_(&ii, &vaux[1]);
4354 if (zu > epsapr[nd] && zv > epsapr[nd]) {
4366 /* ******************************************************************
4368 /* --- OK, the square is valid, the coeff of zone 1 are calculated
4371 minu = (*iordru + 1) << 1;
4373 minv = (*iordrv + 1) << 1;
4376 /* --> If zone 1 is empty, pass to the calculation of Max and Average error. */
4377 if (minu > maxu || minv > maxv) {
4381 /* ----------- The terms connected to degree by U are already calculated -----
4383 /* ------------------- Calculate the coefficients of PATJAC ------------
4386 igsu = minu - i2rdu;
4388 for (jj = minv; jj <= i__1; ++jj) {
4391 for (nd = 1; nd <= i__2; ++nd) {
4392 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4393 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4394 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4395 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4396 patjac_dim2) * patjac_dim1]);
4402 /* --------------- Now the degree is maximally lowered --------
4407 i__1 = 1, i__2 = (*iordru << 1) + 1, i__1 = advapp_max(i__1,i__2);
4408 minu = advapp_max(i__1,*ndminu);
4411 i__1 = 1, i__2 = (*iordrv << 1) + 1, i__1 = advapp_max(i__1,i__2);
4412 minv = advapp_max(i__1,*ndminv);
4416 for (nd = 1; nd <= i__1; ++nd) {
4418 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4419 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4420 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4421 patjac_dim2 * patjac_dim1], &epsapr[nd], &vecerr[
4422 vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4430 /* --> Calculate the average error. */
4431 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4432 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4435 /* --> Set to 0.D0 the rejected coeffs. */
4436 i__2 = idim + ndses - 1;
4437 for (ii = idim; ii <= i__2; ++ii) {
4439 for (jj = nv1; jj <= i__3; ++jj) {
4441 for (kk = nu1; kk <= i__4; ++kk) {
4442 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4451 /* --> Return the nb of coeffs of approximation. */
4452 *ndegpu = advapp_max(*ndegpu,nu);
4453 *ndegpv = advapp_max(*ndegpv,nv);
4458 /* ******************************************************************
4460 /* -------------------- IT IS NOT POSSIBLE TO CUT -------------------
4462 /* ******************************************************************
4466 minu = (*iordru + 1) << 1;
4468 minv = (*iordrv + 1) << 1;
4471 /* ---------------- Calculate the terms connected to the degree by U ------------
4475 for (nd = 1; nd <= i__1; ++nd) {
4477 for (kk = minu; kk <= i__2; ++kk) {
4479 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4480 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4481 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4482 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4483 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4484 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4485 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4489 /* ---------------------- Calculate all coefficients -------
4492 igsu = minu - i2rdu;
4494 for (jj = minv; jj <= i__2; ++jj) {
4496 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4497 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4498 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4499 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4500 patjac_dim2) * patjac_dim1]);
4506 /* ----- Contribution of calculated terms to the approximation error */
4507 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4511 for (nd = 1; nd <= i__1; ++nd) {
4513 minu = (*iordru + 1) << 1;
4517 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4518 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4519 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4523 minv = (*iordrv + 1) << 1;
4526 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4527 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4528 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4532 /* ---------------------------- IF ERRMAX > EPSAPR, stop --------
4535 if (errmax[nd] > epsapr[nd]) {
4540 /* ------------- Otherwise, try to remove again the coeff
4545 i__2 = 1, i__3 = (*iordru << 1) + 1, i__2 = advapp_max(i__2,i__3);
4546 minu = advapp_max(i__2,*ndminu);
4549 i__2 = 1, i__3 = (*iordrv << 1) + 1, i__2 = advapp_max(i__2,i__3);
4550 minv = advapp_max(i__2,*ndminv);
4552 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4553 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &
4554 maxv, iordru, iordrv, xmaxju, xmaxjv, &patjac[
4555 idim * patjac_dim2 * patjac_dim1], &epsapr[nd], &
4556 vecerr[vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4563 /* --------------------- Calculate the average error -------------
4568 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4569 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4572 /* --------------------- Set to 0.D0 the rejected coeffs ----------
4575 i__2 = idim + ndses - 1;
4576 for (ii = idim; ii <= i__2; ++ii) {
4578 for (jj = nv1; jj <= i__3; ++jj) {
4580 for (kk = nu1; kk <= i__4; ++kk) {
4581 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4590 /* --------------- Return the nb of coeff of approximation ---
4593 *ndegpu = advapp_max(*ndegpu,nu);
4594 *ndegpv = advapp_max(*ndegpv,nv);
4602 /* ------------------------------ The end -------------------------------
4604 /* --> Error in inputs */
4609 /* --------- Management of cuts, it is required 0 < NUMDEC <= 5 -------
4612 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4614 if (*numdec <= 0 || *numdec > 5) {
4623 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4625 if (*numdec <= 0 || *numdec > 5) {
4634 /* --> Here it is possible and necessary to cut, choose by 4 if it is possible */
4636 if (*numdec <= 0 || *numdec > 5) {
4641 } else if (*numdec == 2 || *numdec == 4) {
4643 } else if (*numdec == 1 || *numdec == 3) {
4651 AdvApp2Var_SysBase::maermsg_("MMA2CE2", iercod, 7L);
4653 AdvApp2Var_SysBase::mgsomsg_("MMA2CE2", 7L);
4658 //=======================================================================
4659 //function : mma2cfu_
4661 //=======================================================================
4662 int mma2cfu_(integer *ndujac,
4674 /* System generated locals */
4675 integer sosotb_dim1, disotb_dim1, disotb_offset, soditb_dim1,
4676 soditb_offset, diditb_dim1, i__1, i__2;
4678 /* Local variables */
4680 integer nptu2, nptv2, ii, jj;
4681 doublereal bid0, bid1, bid2;
4683 /* **********************************************************************
4688 /* Calculate the terms connected to degree NDUJAC by U of the polynomial approximation */
4689 /* of function F(u,v), starting from its discretisation */
4690 /* on the roots of Legendre polynom of degree */
4691 /* NBPNTU by U and NBPNTV by V. */
4695 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4697 /* INPUT ARGUMENTSE : */
4698 /* ------------------ */
4699 /* NDUJAC: Fixed degree by U for which the terms */
4700 /* allowing to obtain the Legendre or Jacobi coeff*/
4701 /* of even or uneven degree by V are calculated. */
4702 /* NBPNTU: Degree of Legendre polynom on the roots which of */
4703 /* the coefficients of integration by U are calculated */
4704 /* by Gauss method. It is required that NBPNTU = 30, 40, 50 or 61. */
4705 /* NBPNTV: Degree of Legendre polynom on the roots which of */
4706 /* the coefficients of integration by V are calculated */
4707 /* by Gauss method. It is required that NBPNTV = 30, 40, 50 or 61. */
4708 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
4709 /* with ui and vj positive roots of Legendre polynom */
4710 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4711 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
4712 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
4713 /* SOSOTB(0,0) contains F(0,0). */
4714 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
4715 /* with ui and vj positive roots of Legendre polynom */
4716 /* of degree NBPNTU and NBPNTV respectively. */
4717 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
4718 /* with ui and vj positive roots of Legendre polynom */
4719 /* of degree NBPNTU and NBPNTV respectively. */
4720 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
4721 /* avec ui and vj positive roots of Legendre polynom */
4722 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4723 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
4724 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
4725 /* GSSUTB: Table of coefficients of integration by Gauss method */
4726 /* Gauss by U for fixed NDUJAC : i varies from 0 to NBPNTU/2. */
4728 /* OUTPUT ARGUMENTS : */
4729 /* ------------------- */
4730 /* CHPAIR: Table of terms connected to degree NDUJAC by U to calculate the */
4731 /* coeff. of the approximation of EVEN degree by V. */
4732 /* CHIMPR: Table of terms connected to degree NDUJAC by U to calculate */
4733 /* the coeff. of approximation of UNEVEN degree by V. */
4735 /* COMMONS USED : */
4736 /* ---------------- */
4738 /* REFERENCES CALLED : */
4739 /* ----------------------- */
4741 /* DESCRIPTION/NOTES/LIMITATIONS : */
4742 /* ----------------------------------- */
4746 /* **********************************************************************
4748 /* Name of the routine */
4751 /* --------------------------- Initialisations --------------------------
4754 /* Parameter adjustments */
4756 diditb_dim1 = *nbpntu / 2 + 1;
4757 soditb_dim1 = *nbpntu / 2;
4758 soditb_offset = soditb_dim1 + 1;
4759 soditb -= soditb_offset;
4760 disotb_dim1 = *nbpntu / 2;
4761 disotb_offset = disotb_dim1 + 1;
4762 disotb -= disotb_offset;
4763 sosotb_dim1 = *nbpntu / 2 + 1;
4766 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4768 AdvApp2Var_SysBase::mgenmsg_("MMA2CFU", 7L);
4771 nptu2 = *nbpntu / 2;
4772 nptv2 = *nbpntv / 2;
4774 /* **********************************************************************
4776 /* CALCULATE COEFFICIENTS BY U */
4778 /* ----------------- Calculate coefficients of even degree --------------
4781 if (*ndujac % 2 == 0) {
4783 for (jj = 1; jj <= i__1; ++jj) {
4787 for (ii = 1; ii <= i__2; ++ii) {
4789 bid1 += sosotb[ii + jj * sosotb_dim1] * bid0;
4790 bid2 += soditb[ii + jj * soditb_dim1] * bid0;
4798 /* --------------- Calculate coefficients of uneven degree ----------
4803 for (jj = 1; jj <= i__1; ++jj) {
4807 for (ii = 1; ii <= i__2; ++ii) {
4809 bid1 += disotb[ii + jj * disotb_dim1] * bid0;
4810 bid2 += diditb[ii + jj * diditb_dim1] * bid0;
4819 /* ------- Add terms connected to the supplementary root (0.D0) ------ */
4820 /* ----------- of Legendre polynom of uneven degree NBPNTU -----------
4822 /* --> Only even NDUJAC terms are modified as GSSUTB(0) = 0 */
4823 /* when NDUJAC is uneven. */
4825 if (*nbpntu % 2 != 0 && *ndujac % 2 == 0) {
4828 for (jj = 1; jj <= i__1; ++jj) {
4829 chpair[jj] += sosotb[jj * sosotb_dim1] * bid0;
4830 chimpr[jj] += diditb[jj * diditb_dim1] * bid0;
4835 /* ------ Calculate the terms connected to supplementary roots (0.D0) ------
4837 /* ----------- of Legendre polynom of uneven degree NBPNTV -----------
4840 if (*nbpntv % 2 != 0) {
4841 /* --> Only CHPAIR terms are calculated as GSSVTB(0,IH-IDEBV)=0
4843 /* when IH is uneven (see MMA2CFV). */
4845 if (*ndujac % 2 == 0) {
4848 for (ii = 1; ii <= i__1; ++ii) {
4849 bid1 += sosotb[ii] * gssutb[ii];
4856 for (ii = 1; ii <= i__1; ++ii) {
4857 bid1 += diditb[ii] * gssutb[ii];
4862 if (*nbpntu % 2 != 0) {
4863 chpair[0] += sosotb[0] * gssutb[0];
4867 /* ------------------------------ The end -------------------------------
4871 AdvApp2Var_SysBase::mgsomsg_("MMA2CFU", 7L);
4876 //=======================================================================
4877 //function : mma2cfv_
4879 //=======================================================================
4880 int mma2cfv_(integer *ndvjac,
4890 /* System generated locals */
4891 integer chpair_dim1, chpair_offset, chimpr_dim1, chimpr_offset,
4892 patjac_offset, i__1, i__2;
4894 /* Local variables */
4896 integer nptv2, ii, jj;
4899 /* **********************************************************************
4904 /* Calculate the coefficients of polynomial approximation of F(u,v) */
4905 /* of degree NDVJAC by V and of degree by U varying from MINDGU to MAXDGU.
4910 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4912 /* INPUT ARGUMENTS : */
4913 /* ------------------ */
4915 /* NDVJAC: Degree of the polynom of approximation by V. */
4916 /* The representation in the orthogonal base starts from degre 0. */
4917 /* The polynomial base is the base of Jacobi of order -1 */
4918 /* (Legendre), 0, 1 or 2 */
4919 /* MINDGU: Degree minimum by U of coeff. to calculate. */
4920 /* MAXDGU: Degree maximum by U of coeff. to calculate. */
4921 /* NBPNTV: Degree of the Legendre polynom on the roots which of */
4922 /* the coefficients of integration by V are calculated */
4923 /* by Gauss method. It is reqired that NBPNTV = 30, 40, 50 or 61 and NDVJAC < NBPNTV. */
4924 /* GSSVTB: Table of coefficients of integration by Gauss method */
4925 /* by V for NDVJAC fixed: j varies from 0 to NBPNTV/2. */
4926 /* CHPAIR: Table of terms connected to degrees from MINDGU to MAXDGU by U to */
4927 /* calculate the coeff. of approximation of EVEN degree NDVJAC by V. */
4928 /* CHIMPR: Table of terms connected to degrees from MINDGU to MAXDGU by U to */
4929 /* calculate the coeff. of approximation of UNEVEN degree NDVJAC by V. */
4931 /* OUTPUT ARGUMENTS : */
4932 /* ------------------- */
4933 /* PATJAC: Table of coefficients by U of the polynom of approximation */
4934 /* P(u,v) of degree MINDGU to MAXDGU by U and NDVJAC by V. */
4936 /* COMMONS USED : */
4937 /* -------------- */
4939 /* REFERENCES CALLED : */
4940 /* --------------------- */
4942 /* DESCRIPTION/NOTES/LIMITATIONS : */
4943 /* ------------------------------- */
4945 /* **********************************************************************
4947 /* Name of the routine */
4950 /* --------------------------- Initialisations --------------------------
4953 /* Parameter adjustments */
4954 patjac_offset = *mindgu;
4955 patjac -= patjac_offset;
4956 chimpr_dim1 = *nbpntv / 2;
4957 chimpr_offset = chimpr_dim1 * *mindgu + 1;
4958 chimpr -= chimpr_offset;
4959 chpair_dim1 = *nbpntv / 2 + 1;
4960 chpair_offset = chpair_dim1 * *mindgu;
4961 chpair -= chpair_offset;
4964 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4966 AdvApp2Var_SysBase::mgenmsg_("MMA2CFV", 7L);
4968 nptv2 = *nbpntv / 2;
4970 /* --------- Calculate the coefficients for even degree NDVJAC ----------
4973 if (*ndvjac % 2 == 0) {
4975 for (ii = *mindgu; ii <= i__1; ++ii) {
4978 for (jj = 1; jj <= i__2; ++jj) {
4979 bid1 += chpair[jj + ii * chpair_dim1] * gssvtb[jj];
4986 /* -------- Calculate the coefficients for uneven degree NDVJAC -----
4991 for (ii = *mindgu; ii <= i__1; ++ii) {
4994 for (jj = 1; jj <= i__2; ++jj) {
4995 bid1 += chimpr[jj + ii * chimpr_dim1] * gssvtb[jj];
5003 /* ------- Add terms connected to the supplementary root (0.D0) ----- */
5004 /* --------of the Legendre polynom of uneven degree NBPNTV --------- */
5006 if (*nbpntv % 2 != 0 && *ndvjac % 2 == 0) {
5009 for (ii = *mindgu; ii <= i__1; ++ii) {
5010 patjac[ii] += bid1 * chpair[ii * chpair_dim1];
5015 /* ------------------------------ The end -------------------------------
5019 AdvApp2Var_SysBase::mgsomsg_("MMA2CFV", 7L);
5024 //=======================================================================
5025 //function : mma2ds1_
5027 //=======================================================================
5028 int AdvApp2Var_ApproxF2var::mma2ds1_(integer *ndimen,
5031 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
5046 /* System generated locals */
5047 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5048 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5049 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5050 fpntab_offset, i__1;
5052 /* Local variables */
5054 integer ibid1, ibid2, iuouv, nd;
5057 /* **********************************************************************
5062 /* Discretisation of function F(u,v) on the roots of Legendre polynoms. */
5066 /* FONCTION&,DISCRETISATION,&POINT */
5068 /* INPUT ARGUMENTS : */
5069 /* ------------------ */
5070 /* NDIMEN: Dimension of the space. */
5071 /* UINTFN: Limits of the interval of definition by u of the function */
5072 /* to be processed: (UINTFN(1),UINTFN(2)). */
5073 /* VINTFN: Limits of the interval of definition by v of the function */
5074 /* to be processed: (VINTFN(1),VINTFN(2)). */
5075 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5076 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5077 /* FONCNP is discretized by u. */
5078 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5079 /* FONCNP is discretized by v. */
5080 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5081 /* of Legendre of degree NBPNTU defined on (-1,1). */
5082 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5083 /* of Legendre of degree NBPNTV defined on (-1,1). */
5084 /* ISOFAV: Shows the type of iso of F(u,v) to be extracted to improve */
5085 /* the rapidity of calculation (has no influence on the form */
5087 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5088 /* with fixed u (with NBPNTV values different from v). */
5089 /* = 2, shows that it is necessaty to calculate the points of F(u,v) */
5090 /* with fixed v (with NBPNTU values different from u). */
5091 /* SOSOTB: Preinitialized table (input/output argument). */
5092 /* DISOTB: Preinitialized table (input/output argument). */
5093 /* SODITB: Preinitialized table (input/output argument). */
5094 /* DIDITB: Preinitialized table (input/output argument). */
5096 /* OUTPUT ARGUMENTS : */
5097 /* ------------------- */
5098 /* SOSOTB: Table where the terms */
5099 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5100 /* are added with ui and vj positive roots of Legendre polynom */
5101 /* of degree NBPNTU and NBPNTV respectively. */
5102 /* DISOTB: Table where the terms */
5103 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5104 /* are added with ui and vj positive roots of Legendre polynom */
5105 /* of degree NBPNTU and NBPNTV respectively. */
5106 /* SODITB: Table where the terms */
5107 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5108 /* are added with ui and vj positive roots of Legendre polynom */
5109 /* of degree NBPNTU and NBPNTV respectively. */
5110 /* DIDITB: Table where the terms */
5111 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5112 /* are added with ui and vj positive roots of Legendre polynom */
5113 /* of degree NBPNTU and NBPNTV respectively. */
5114 /* FPNTAB: Auxiliary table. */
5115 /* TTABLE: Auxiliary table. */
5116 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5117 /* the returned error code is equal to error code of FONCNP + 100. */
5119 /* COMMONS USED : */
5120 /* ---------------- */
5122 /* REFERENCES CALLED : */
5123 /* --------------------- */
5125 /* DESCRIPTION/NOTES/LIMITATIONS : */
5126 /* ----------------------------------- */
5127 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5128 /* where MA2FXK should be in the following form : */
5129 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,ISOFAV,TCONST,NBPTAB */
5130 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5131 /* with the following input arguments : */
5132 /* - NDIMEN is integer defined as the sum of dimensions of */
5133 /* sub-spaces (i.e. total dimension of the problem). */
5134 /* - UINTFN(2) is a table of 2 reals containing the interval */
5135 /* by u where the function to be approximated is defined */
5136 /* (so it is equal to UIFONC). */
5137 /* - VINTFN(2) is a table of 2 reals containing the interval */
5138 /* by v where the function to be approximated is defined */
5139 /* (so it is equal to VIFONC). */
5140 /* - ISOFAV, is 1 if it is necessary to calculate points with constant u, */
5141 /* is 2 if it is necessary to calculate points with constant v. */
5142 /* Any other value is an error. */
5143 /* - TCONST, real, value of the fixed parameter. Takes values */
5144 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5145 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5146 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5147 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5148 /* 'free' parameter of discretization (v if IISOFAV=1, */
5149 /* u if IISOFAV=2). */
5150 /* - IDERIU, integer, takes values between 0 (position) */
5151 /* and IORDRE(1) (partial derivative of the function by u */
5152 /* of order IORDRE(1) if IORDRE(1) > 0). */
5153 /* - IDERIV, integer, takes values between 0 (position) */
5154 /* and IORDRE(2) (partial derivative of the function by v */
5155 /* of order IORDRE(2) if IORDRE(2) > 0). */
5156 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5157 /* points of the derivative : */
5164 /* and the output arguments aret : */
5165 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5166 /* NBPTAB points calculated in FONCNP. */
5167 /* - IERCOD is, at output the error code of FONCNP. This code */
5168 /* (integer) should be strictly positive if there is a problem. */
5170 /* The input arguments SHOULD NOT be modified under FONCNP.
5173 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5174 /* values of UROOTB and VROOTB are consequently modified. */
5176 /* -->The results of discretisation are ranked in 4 tables */
5177 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5178 /* during the calculation of coefficients of the polynom of approximation. */
5180 /* When NBPNTU is uneven : */
5181 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5182 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5183 /* When NBPNTV is uneven : */
5184 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5185 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5186 /* When NBPNTU and NBPNTV are uneven : */
5187 /* term SOSOTB(0,0) contains F(0,0). */
5190 /* **********************************************************************
5192 /* Name of the routine */
5195 /* --------------------------- Initialization --------------------------
5198 /* Parameter adjustments */
5199 fpntab_dim1 = *ndimen;
5200 fpntab_offset = fpntab_dim1 + 1;
5201 fpntab -= fpntab_offset;
5205 diditb_dim1 = *nbpntu / 2 + 1;
5206 diditb_dim2 = *nbpntv / 2 + 1;
5207 diditb_offset = diditb_dim1 * diditb_dim2;
5208 diditb -= diditb_offset;
5209 soditb_dim1 = *nbpntu / 2;
5210 soditb_dim2 = *nbpntv / 2;
5211 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5212 soditb -= soditb_offset;
5213 disotb_dim1 = *nbpntu / 2;
5214 disotb_dim2 = *nbpntv / 2;
5215 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5216 disotb -= disotb_offset;
5217 sosotb_dim1 = *nbpntu / 2 + 1;
5218 sosotb_dim2 = *nbpntv / 2 + 1;
5219 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5220 sosotb -= sosotb_offset;
5225 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5227 AdvApp2Var_SysBase::mgenmsg_("MMA2DS1", 7L);
5230 if (*isofav < 1 || *isofav > 2) {
5236 /* **********************************************************************
5238 /* --------- Discretization by U on the roots of the polynom of ------ */
5239 /* --------------- Legendre of degree NBPNTU, iso-V by iso-V --------- */
5240 /* **********************************************************************
5244 mma2ds2_(ndimen, &uintfn[1], &vintfn[1], foncnp, nbpntu, nbpntv, &
5245 urootb[1], &vrootb[1], &iuouv, &sosotb[sosotb_offset], &
5246 disotb[disotb_offset], &soditb[soditb_offset], &diditb[
5247 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5249 /* ******************************************************************
5251 /* --------- Discretization by V on the roots of the polynom of ------ */
5252 /* --------------- Legendre of degree NBPNTV, iso-V by iso-V --------- */
5253 /* ******************************************************************
5257 /* --> Inversion of indices of tables */
5259 for (nd = 1; nd <= i__1; ++nd) {
5260 isz1 = *nbpntu / 2 + 1;
5261 isz2 = *nbpntv / 2 + 1;
5262 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5263 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5264 ibid1, &ibid2, iercod);
5268 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5269 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5270 ibid1, &ibid2, iercod);
5276 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5277 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5278 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5282 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5283 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5284 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5291 mma2ds2_(ndimen, &vintfn[1], &uintfn[1], foncnp, nbpntv, nbpntu, &
5292 vrootb[1], &urootb[1], &iuouv, &sosotb[sosotb_offset], &
5293 soditb[soditb_offset], &disotb[disotb_offset], &diditb[
5294 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5295 /* --> Inversion of indices of tables */
5297 for (nd = 1; nd <= i__1; ++nd) {
5298 isz1 = *nbpntv / 2 + 1;
5299 isz2 = *nbpntu / 2 + 1;
5300 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5301 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5302 ibid1, &ibid2, iercod);
5306 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5307 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5308 ibid1, &ibid2, iercod);
5314 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5315 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5316 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5320 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5321 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5322 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5330 /* ------------------------------ The end -------------------------------
5336 AdvApp2Var_SysBase::maermsg_("MMA2DS1", iercod, 7L);
5339 AdvApp2Var_SysBase::mgsomsg_("MMA2DS1", 7L);
5344 //=======================================================================
5345 //function : mma2ds2_
5347 //=======================================================================
5348 int mma2ds2_(integer *ndimen,
5351 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
5367 /* System generated locals */
5368 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5369 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5370 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5371 fpntab_offset, i__1, i__2, i__3;
5373 /* Local variables */
5376 doublereal alinu, blinu, alinv, blinv, tcons;
5377 doublereal dbfn1[2], dbfn2[2];
5378 integer nuroo, nvroo, id, iu, iv;
5382 /* **********************************************************************
5387 /* Discretization of function F(u,v) on the roots of polynoms of Legendre. */
5391 /* FONCTION&,DISCRETISATION,&POINT */
5393 /* INPUT ARGUMENTS : */
5394 /* ------------------ */
5395 /* NDIMEN: Dimension of the space. */
5396 /* UINTFN: Limits of the interval of definition by u of the function */
5397 /* to be processed: (UINTFN(1),UINTFN(2)). */
5398 /* VINTFN: Limits of the interval of definition by v of the function */
5399 /* to be processed: (VINTFN(1),VINTFN(2)). */
5400 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5401 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5402 /* FONCNP is discretized by u. */
5403 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5404 /* FONCNP is discretized by v. */
5405 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5406 /* of Legendre of degree NBPNTU defined on (-1,1). */
5407 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5408 /* of Legendre of degree NBPNTV defined on (-1,1). */
5409 /* IIUOUV: Shows the type of iso of F(u,v) tom be extracted to improve the */
5410 /* rapidity of calculation (has no influence on the form of result) */
5411 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5412 /* with fixed u (so with NBPNTV values different from v). */
5413 /* = 2, shows that it is necessary to calculate the points of F(u,v) */
5414 /* with fixed v (so with NBPNTV values different from u). */
5415 /* SOSOTB: Preinitialized table (input/output argument). */
5416 /* DISOTB: Preinitialized table (input/output argument). */
5417 /* SODITB: Preinitialized table (input/output argument). */
5418 /* DIDITB: Preinitialized table (input/output argument). */
5420 /* OUTPUT ARGUMENTS : */
5421 /* ------------------- */
5422 /* SOSOTB: Table where the terms */
5423 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5424 /* are added with ui and vj positive roots of Legendre polynom */
5425 /* of degree NBPNTU and NBPNTV respectively. */
5426 /* DISOTB: Table where the terms */
5427 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5428 /* are added with ui and vj positive roots of Legendre polynom */
5429 /* of degree NBPNTU and NBPNTV respectively. */
5430 /* SODITB: Table where the terms */
5431 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5432 /* are added with ui and vj positive roots of Legendre polynom */
5433 /* of degree NBPNTU and NBPNTV respectively. */
5434 /* DIDITB: Table where the terms */
5435 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5436 /* are added with ui and vj positive roots of Legendre polynom */
5437 /* of degree NBPNTU and NBPNTV respectively. */
5438 /* FPNTAB: Auxiliary table. */
5439 /* TTABLE: Auxiliary table. */
5440 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5441 /* the returned error code is equal to error code of FONCNP + 100. */
5443 /* COMMONS USED : */
5444 /* ---------------- */
5446 /* REFERENCES CALLED : */
5447 /* --------------------- */
5449 /* DESCRIPTION/NOTES/LIMITATIONS : */
5450 /* ----------------------------------- */
5451 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5452 /* where MA2FXK should be in the following form : */
5453 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIIUOUV,TCONST,NBPTAB */
5454 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5455 /* with the following input arguments : */
5456 /* - NDIMEN is integer defined as the sum of dimensions of */
5457 /* sub-spaces (i.e. total dimension of the problem). */
5458 /* - UINTFN(2) is a table of 2 reals containing the interval */
5459 /* by u where the function to be approximated is defined */
5460 /* (so it is equal to UIFONC). */
5461 /* - VINTFN(2) is a table of 2 reals containing the interval */
5462 /* by v where the function to be approximated is defined */
5463 /* (so it is equal to VIFONC). */
5464 /* - IIIUOUV, is 1 if it is necessary to calculate points with constant u, */
5465 /* is 2 if it is necessary to calculate points with constant v. */
5466 /* Any other value is an error. */
5467 /* - TCONST, real, value of the fixed parameter. Takes values */
5468 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5469 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5470 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5471 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5472 /* 'free' parameter of discretization (v if IIIUOUV=1, */
5473 /* u if IIIUOUV=2). */
5474 /* - IDERIU, integer, takes values between 0 (position) */
5475 /* and IORDRE(1) (partial derivative of the function by u */
5476 /* of order IORDRE(1) if IORDRE(1) > 0). */
5477 /* - IDERIV, integer, takes values between 0 (position) */
5478 /* and IORDRE(2) (partial derivative of the function by v */
5479 /* of order IORDRE(2) if IORDRE(2) > 0). */
5480 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5481 /* points of the derivative : */
5488 /* and the output arguments aret : */
5489 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5490 /* NBPTAB points calculated in FONCNP. */
5491 /* - IERCOD is, at output the error code of FONCNP. This code */
5492 /* (integer) should be strictly positive if there is a problem. */
5494 /* The input arguments SHOULD NOT be modified under FONCNP.
5497 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5498 /* values of UROOTB and VROOTB are consequently modified. */
5500 /* -->The results of discretisation are ranked in 4 tables */
5501 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5502 /* during the calculation of coefficients of the polynom of approximation. */
5504 /* When NBPNTU is uneven : */
5505 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5506 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5507 /* When NBPNTV is uneven : */
5508 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5509 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5510 /* When NBPNTU and NBPNTV are uneven : */
5511 /* term SOSOTB(0,0) contains F(0,0). */
5513 /* ATTENTION: These 4 tables are filled by varying the */
5514 /* 1st index first. So, the discretizations */
5515 /* of F(...,t) (for IIUOUV = 2) or of F(t,...) (IIUOUV = 1) */
5516 /* are stored in SOSOTB(...,t), SODITB(...,t), etc... */
5517 /* (this allows to gain important time). */
5518 /* It is required that the caller, in case of IIUOUV=1, */
5519 /* invert the roles of u and v, of SODITB and DISOTB BEFORE the */
5522 /* **********************************************************************
5525 /* Name of the routine */
5527 /* --> Indices of loops. */
5529 /* --------------------------- Initialization --------------------------
5532 /* Parameter adjustments */
5536 fpntab_dim1 = *ndimen;
5537 fpntab_offset = fpntab_dim1 + 1;
5538 fpntab -= fpntab_offset;
5540 diditb_dim1 = *nbpntu / 2 + 1;
5541 diditb_dim2 = *nbpntv / 2 + 1;
5542 diditb_offset = diditb_dim1 * diditb_dim2;
5543 diditb -= diditb_offset;
5544 soditb_dim1 = *nbpntu / 2;
5545 soditb_dim2 = *nbpntv / 2;
5546 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5547 soditb -= soditb_offset;
5548 disotb_dim1 = *nbpntu / 2;
5549 disotb_dim2 = *nbpntv / 2;
5550 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5551 disotb -= disotb_offset;
5552 sosotb_dim1 = *nbpntu / 2 + 1;
5553 sosotb_dim2 = *nbpntv / 2 + 1;
5554 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5555 sosotb -= sosotb_offset;
5559 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5561 AdvApp2Var_SysBase::mgenmsg_("MMA2DS2", 7L);
5565 alinu = (uintfn[2] - uintfn[1]) / 2.;
5566 blinu = (uintfn[2] + uintfn[1]) / 2.;
5567 alinv = (vintfn[2] - vintfn[1]) / 2.;
5568 blinv = (vintfn[2] + vintfn[1]) / 2.;
5571 dbfn1[0] = vintfn[1];
5572 dbfn1[1] = vintfn[2];
5573 dbfn2[0] = uintfn[1];
5574 dbfn2[1] = uintfn[2];
5576 dbfn1[0] = uintfn[1];
5577 dbfn1[1] = uintfn[2];
5578 dbfn2[0] = vintfn[1];
5579 dbfn2[1] = vintfn[2];
5582 /* **********************************************************************
5584 /* -------- Discretization by U on the roots of Legendre polynom -------- */
5585 /* ---------------- of degree NBPNTU, with Vj fixed -------------------- */
5586 /* **********************************************************************
5589 nuroo = *nbpntu / 2;
5590 nvroo = *nbpntv / 2;
5591 jdec = (*nbpntu + 1) / 2;
5593 /* ----------- Loading of parameters of discretization by U ------------- */
5596 for (iu = 1; iu <= i__1; ++iu) {
5597 ttable[iu] = blinu + alinu * urootb[iu];
5601 /* -------------- For Vj fixed, negative root of Legendre ------------- */
5604 for (iv = 1; iv <= i__1; ++iv) {
5605 tcons = blinv + alinv * vrootb[iv];
5606 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5607 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5608 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5613 for (id = 1; id <= i__2; ++id) {
5615 for (iu = 1; iu <= i__3; ++iu) {
5616 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5617 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5618 sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1]
5619 = sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) *
5620 sosotb_dim1] + up + um;
5621 disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) * disotb_dim1]
5622 = disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) *
5623 disotb_dim1] + up - um;
5624 soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) * soditb_dim1]
5625 = soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) *
5626 soditb_dim1] - up - um;
5627 diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1]
5628 = diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) *
5629 diditb_dim1] - up + um;
5632 if (*nbpntu % 2 != 0) {
5633 up = fpntab[id + jdec * fpntab_dim1];
5634 sosotb[(nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1] +=
5636 diditb[(nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1] -=
5644 /* --------- For Vj = 0 (uneven NBPNTV), discretization by U ----------- */
5646 if (*nbpntv % 2 != 0) {
5648 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5649 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5650 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5655 for (id = 1; id <= i__1; ++id) {
5657 for (iu = 1; iu <= i__2; ++iu) {
5658 up = fpntab[id + (jdec + iu) * fpntab_dim1];
5659 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5660 sosotb[iu + id * sosotb_dim2 * sosotb_dim1] = sosotb[iu + id *
5661 sosotb_dim2 * sosotb_dim1] + up + um;
5662 diditb[iu + id * diditb_dim2 * diditb_dim1] = diditb[iu + id *
5663 diditb_dim2 * diditb_dim1] + up - um;
5666 if (*nbpntu % 2 != 0) {
5667 up = fpntab[id + jdec * fpntab_dim1];
5668 sosotb[id * sosotb_dim2 * sosotb_dim1] += up;
5674 /* -------------- For Vj fixed, positive root of Legendre ------------- */
5677 for (iv = 1; iv <= i__1; ++iv) {
5678 tcons = alinv * vrootb[(*nbpntv + 1) / 2 + iv] + blinv;
5679 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5680 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5681 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5686 for (id = 1; id <= i__2; ++id) {
5688 for (iu = 1; iu <= i__3; ++iu) {
5689 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5690 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5691 sosotb[iu + (iv + id * sosotb_dim2) * sosotb_dim1] = sosotb[
5692 iu + (iv + id * sosotb_dim2) * sosotb_dim1] + up + um;
5693 disotb[iu + (iv + id * disotb_dim2) * disotb_dim1] = disotb[
5694 iu + (iv + id * disotb_dim2) * disotb_dim1] + up - um;
5695 soditb[iu + (iv + id * soditb_dim2) * soditb_dim1] = soditb[
5696 iu + (iv + id * soditb_dim2) * soditb_dim1] + up + um;
5697 diditb[iu + (iv + id * diditb_dim2) * diditb_dim1] = diditb[
5698 iu + (iv + id * diditb_dim2) * diditb_dim1] + up - um;
5701 if (*nbpntu % 2 != 0) {
5702 up = fpntab[id + jdec * fpntab_dim1];
5703 sosotb[(iv + id * sosotb_dim2) * sosotb_dim1] += up;
5704 diditb[(iv + id * diditb_dim2) * diditb_dim1] += up;
5711 /* ------------------------------ The end -------------------------------
5717 AdvApp2Var_SysBase::maermsg_("MMA2DS2", iercod, 7L);
5720 AdvApp2Var_SysBase::mgsomsg_("MMA2DS2", 7L);
5725 //=======================================================================
5726 //function : mma2er1_
5728 //=======================================================================
5729 int mma2er1_(integer *ndjacu,
5745 /* System generated locals */
5746 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
5749 /* Local variables */
5754 doublereal bid0, bid1;
5756 /* **********************************************************************
5761 /* Calculate max approximation error done when */
5762 /* the coefficients of PATJAC such that the degree by U varies between */
5763 /* MINDGU and MAXDGU and the degree by V varies between MINDGV and MAXDGV are removed. */
5767 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5769 /* INPUT ARGUMENTS : */
5770 /* ------------------ */
5771 /* NDJACU: Dimension by U of table PATJAC. */
5772 /* NDJACV: Dimension by V of table PATJAC. */
5773 /* NDIMEN: Dimension of the space. */
5774 /* MINDGU: Lower limit of index by U of coeff. of PATJAC to be taken into account. */
5775 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5776 /* MINDGV: Lower limit of index by V of coeff. of PATJAC to be taken into account. */
5777 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5778 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5779 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5780 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5781 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5782 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5783 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5784 /* PATJAC: Table of coeff. of square of approximation with */
5785 /* constraints of order IORDRU by U and IORDRV by V. */
5786 /* VECERR: Auxiliary vector. */
5787 /* ERREUR: MAX Error commited during removal of ALREADY CALCULATED coeff of PATJAC */
5789 /* OUTPUT ARGUMENTS : */
5790 /* ------------------- */
5791 /* ERREUR: MAX Error commited during removal of coeff of PATJAC */
5792 /* of indices from MINDGU to MAXDGU by U and from MINDGV to MAXDGV by V */
5793 /* THEN the already calculated error. */
5795 /* COMMONS USED : */
5796 /* ---------------- */
5798 /* REFERENCES CALLED : */
5799 /* --------------------- */
5801 /* DESCRIPTION/NOTES/LIMITATIONS : */
5802 /* ----------------------------------- */
5803 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5804 /* approximation of F(U,V). The indices i and j show the degree */
5805 /* by U and by V of base polynoms. These polynoms have the form: */
5807 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5809 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5810 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5812 /* The contribution to the error of term Cij when it is */
5813 /* removed from PATJAC is increased by: */
5815 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5817 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5819 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5823 /* ***********************************************************************
5825 /* Name of the routine */
5828 /* ----------------------------- Initialisations ------------------------
5831 /* Parameter adjustments */
5833 patjac_dim1 = *ndjacu + 1;
5834 patjac_dim2 = *ndjacv + 1;
5835 patjac_offset = patjac_dim1 * patjac_dim2;
5836 patjac -= patjac_offset;
5839 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5841 AdvApp2Var_SysBase::mgenmsg_("MMA2ER1", 7L);
5844 minu = (*iordru + 1) << 1;
5845 minv = (*iordrv + 1) << 1;
5847 /* ------------------- Calculate the increment of the max error --------------- */
5848 /* ----- during the removal of the coeffs of indices from MINDGU to MAXDGU ---- */
5849 /* ---------------- by U and indices from MINDGV to MAXDGV by V --------------- */
5852 for (nd = 1; nd <= i__1; ++nd) {
5855 for (jj = *mindgv; jj <= i__2; ++jj) {
5858 for (ii = *mindgu; ii <= i__3; ++ii) {
5859 bid0 += (d__1 = patjac[ii + (jj + nd * patjac_dim2) *
5860 patjac_dim1], advapp_abs(d__1)) * xmaxju[ii - minu];
5863 bid1 = bid0 * xmaxjv[jj - minv] + bid1;
5871 /* ----------------------- Calculate the max error ----------------------*/
5873 bid1 = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
5877 *erreur = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
5879 /* ------------------------- The end ------------------------------------
5883 AdvApp2Var_SysBase::mgsomsg_("MMA2ER1", 7L);
5888 //=======================================================================
5889 //function : mma2er2_
5891 //=======================================================================
5892 int mma2er2_(integer *ndjacu,
5904 doublereal *epmscut,
5911 /* System generated locals */
5912 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2;
5915 /* Local variables */
5918 integer i2rdu, i2rdv;
5919 doublereal errnu, errnv;
5920 integer ii, nd, jj, nu, nv;
5921 doublereal bid0, bid1;
5923 /* **********************************************************************
5928 /* Remove coefficients of PATJAC to obtain the minimum degree */
5929 /* by U and V checking the imposed tolerance. */
5933 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5935 /* INPUT ARGUMENTS : */
5936 /* ------------------ */
5937 /* NDJACU: Degree by U of table PATJAC. */
5938 /* NDJACV: Degree by V of table PATJAC. */
5939 /* NDIMEN: Dimension of the space. */
5940 /* MINDGU: Limit of index by U of coeff. of PATJAC to be PRESERVED (should be >=0). */
5941 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5942 /* MINDGV: Limit of index by V of coeff. of PATJAC to be PRESERVED (should be >=0). */
5943 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5944 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5945 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5946 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5947 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5948 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5949 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5950 /* PATJAC: Table of coeff. of square of approximation with */
5951 /* constraints of order IORDRU by U and IORDRV by V. */
5952 /* EPMSCUT: Tolerance of approximation. */
5953 /* VECERR: Auxiliary vector. */
5954 /* ERREUR: MAX Error commited ALREADY CALCULATED */
5956 /* OUTPUT ARGUMENTS : */
5957 /* ------------------- */
5958 /* ERREUR: MAX Error commited by preserving only coeff of PATJAC */
5959 /* of indices from 0 to NEWDGU by U and from 0 to NEWDGV by V */
5960 /* PLUS the already calculated error. */
5961 /* NEWDGU: Min. Degree by U such as the square of approximation */
5962 /* could check the tolerance. There is always NEWDGU >= MINDGU >= 0. */
5963 /* NEWDGV: Min. Degree by V such as the square of approximation */
5964 /* could check the tolerance. There is always NEWDGV >= MINDGV >= 0. */
5967 /* COMMONS USED : */
5968 /* ---------------- */
5970 /* REFERENCES CALLED : */
5971 /* --------------------- */
5973 /* DESCRIPTION/NOTES/LIMITATIONS : */
5974 /* ----------------------------------- */
5975 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5976 /* approximation of F(U,V). The indices i and j show the degree */
5977 /* by U and by V of base polynoms. These polynoms have the form: */
5979 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5981 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5982 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5984 /* The contribution to the error of term Cij when it is */
5985 /* removed from PATJAC is increased by: */
5987 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5989 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5991 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5995 /* **********************************************************************
5997 /* Name of the routine */
6000 /* ----------------------------- Initialisations ------------------------
6003 /* Parameter adjustments */
6005 patjac_dim1 = *ndjacu + 1;
6006 patjac_dim2 = *ndjacv + 1;
6007 patjac_offset = patjac_dim1 * patjac_dim2;
6008 patjac -= patjac_offset;
6011 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
6013 AdvApp2Var_SysBase::mgenmsg_("MMA2ER2", 7L);
6016 i2rdu = (*iordru + 1) << 1;
6017 i2rdv = (*iordrv + 1) << 1;
6021 /* **********************************************************************
6023 /* -------------------- Cutting of oefficients ------------------------
6025 /* **********************************************************************
6030 /* ------------------- Calculate the increment of max error --------------- */
6031 /* ----- during the removal of coeff. of indices from MINDGU to MAXDGU ------ */
6032 /* ---------------- by U, the degree by V is fixed to NV -----------------
6037 bid0 = xmaxjv[nv - i2rdv];
6039 for (nd = 1; nd <= i__1; ++nd) {
6042 for (ii = i2rdu; ii <= i__2; ++ii) {
6043 bid1 += (d__1 = patjac[ii + (nv + nd * patjac_dim2) *
6044 patjac_dim1], advapp_abs(d__1)) * xmaxju[ii - i2rdu] * bid0;
6051 vecerr[1] = *epmscut * 2;
6053 errnv = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6055 /* ------------------- Calculate the increment of max error --------------- */
6056 /* ----- during the removal of coeff. of indices from MINDGV to MAXDGV ------ */
6057 /* ---------------- by V, the degree by U is fixed to NU -----------------
6062 bid0 = xmaxju[nu - i2rdu];
6064 for (nd = 1; nd <= i__1; ++nd) {
6067 for (jj = i2rdv; jj <= i__2; ++jj) {
6068 bid1 += (d__1 = patjac[nu + (jj + nd * patjac_dim2) *
6069 patjac_dim1], advapp_abs(d__1)) * xmaxjv[jj - i2rdv] * bid0;
6076 vecerr[1] = *epmscut * 2;
6078 errnu = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6080 /* ----------------------- Calculate the max error ----------------------
6086 errnu = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6088 errnv = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6090 if (errnu > errnv) {
6091 if (errnv < *epmscut) {
6098 if (errnu < *epmscut) {
6108 /* -------------------------- Return the degrees -------------------
6112 *newdgu = advapp_max(nu,1);
6113 *newdgv = advapp_max(nv,1);
6115 /* ----------------------------------- The end --------------------------
6119 AdvApp2Var_SysBase::mgsomsg_("MMA2ER2", 7L);
6124 //=======================================================================
6125 //function : mma2fnc_
6127 //=======================================================================
6128 int AdvApp2Var_ApproxF2var::mma2fnc_(integer *ndimen,
6132 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
6158 /* System generated locals */
6159 integer courbe_dim1, courbe_dim2, courbe_offset, somtab_dim1, somtab_dim2,
6160 somtab_offset, diftab_dim1, diftab_dim2, diftab_offset,
6161 contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
6162 contr2_offset, errmax_dim1, errmax_offset, errmoy_dim1,
6163 errmoy_offset, i__1;
6166 /* Local variables */
6169 integer ideb1, ibid1, ibid2, ncfja, ndgre, ilong,
6171 doublereal* wrkar = 0;
6174 doublereal uvpav[4] /* was [2][2] */;
6178 doublereal uv11[4] /* was [2][2] */;
6181 integer isz1, isz2, isz3, isz4, isz5;
6182 intptr_t ipt1, ipt2, ipt3, ipt4, ipt5,iptt, jptt;
6184 /* **********************************************************************
6189 /* Approximation of a limit of non polynomial function F(u,v) */
6190 /* (in the space of dimension NDIMEN) by SEVERAL */
6191 /* polynomial curves, by the method of least squares. The parameter of the function is preserved. */
6195 /* TOUS, AB_SPECIFI :: FONCTION&,EXTREMITE&, APPROXIMATION, &COURBE. */
6197 /* INPUT ARGUMENTS : */
6198 /* ----------------- */
6199 /* NDIMEN: Total Dimension of the space (sum of dimensions */
6200 /* of sub-spaces) */
6201 /* NBSESP: Number of "independent" sub-spaces. */
6202 /* NDIMSE: Table of dimensions of sub-spaces. */
6203 /* UVFONC: Limits of the interval (a,b)x(c,d) of definition of the */
6204 /* function to be approached by U (UVFONC(*,1) contains (a,b)) */
6205 /* and by V (UVFONC(*,2) contains (c,d)). */
6206 /* FONCNP: External function of position on the non polynomial function to be approached. */
6207 /* TCONST: Value of isoparameter of F(u,v) to be discretized. */
6208 /* ISOFAV: Type of chosen iso, = 1, shose that discretization is with u */
6209 /* fixed; = 2, shows that v is fixed. */
6210 /* NBROOT: Nb of points of discretisation of the iso, extremities not included. */
6211 /* ROOTLG: Table of roots of the polynom of Legendre defined on */
6212 /* (-1,1), of degree NBROOT. */
6213 /* IORDRE: Order of constraint at the extremities of the limit */
6214 /* -1 = no constraints, */
6215 /* 0 = constraints of passage to limits (i.e. C0), */
6216 /* 1 = C0 + constraints of 1st derivatives (i.e. C1), */
6217 /* 2 = C1 + constraints of 2nd derivatives (i.e. C2). */
6218 /* IDERIV: Order of derivative of the limit. */
6219 /* NDGJAC: Degree of serial development to be used for calculation in */
6220 /* the Jacobi base. */
6221 /* NBCRMX: Max Nb of curves to be created. */
6222 /* NCFLIM: Max Nb of coeff of the polynomial curve */
6223 /* of approximation (should be above or equal to */
6224 /* 2*IORDRE+2 and below or equal to 50). */
6225 /* EPSAPR: Table of required errors of approximation */
6226 /* sub-space by sub-space. */
6228 /* OUTPUT ARGUMENTS : */
6229 /* ------------------- */
6230 /* NCOEFF: Number of significative coeff of calculated curves. */
6231 /* COURBE: Table of coeff. of calculated polynomial curves. */
6232 /* Should be dimensioned in (NCFLIM,NDIMEN,NBCRMX). */
6233 /* These curves are ALWAYS parametrized in (-1,1). */
6234 /* NBCRBE: Nb of calculated curves. */
6235 /* SOMTAB: For F defined on (-1,1) (otherwise rescale the */
6236 /* parameters), this is the table of sums F(u,vj) + F(u,-vj)
6238 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6239 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6240 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6241 /* v of order IDERIV are taken). */
6242 /* DIFTAB: For F defined on (-1,1) (otherwise rescale the */
6243 /* parameters), this is the table of sums F(u,vj) - F(u,-vj)
6245 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6246 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6247 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6248 /* v of order IDERIV are taken). */
6249 /* CONTR1: Contains the coordinates of the left extremity of the iso */
6250 /* and of its derivatives till order IORDRE */
6251 /* CONTR2: Contains the coordinates of the right extremity of the iso */
6252 /* and of its derivatives till order IORDRE */
6253 /* TABDEC: Table of NBCRBE+1 parameters of cut of UVFONC(1:2,1)
6255 /* if ISOFAV=2, or of UVFONC(1:2,2) if ISOFAV=1. */
6256 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6257 /* committed in the approximation of FONCNP by NBCRBE curves. */
6258 /* ERRMOY: Table of AVERAGE errors (sub-space by sub-space) */
6259 /* committed in the approximation of FONCNP by NBCRBE curves. */
6260 /* IERCOD: Error code: */
6261 /* -1 = ERRMAX > EPSAPR for at least one sub-space. */
6262 /* (the resulting curves of at least mathematic degree NCFLIM-1 */
6263 /* are calculated). */
6264 /* 0 = Everything is ok. */
6265 /* 1 = Pb of incoherence of inputs. */
6266 /* 10 = Pb of calculation of the interpolation of constraints. */
6267 /* 13 = Pb in the dynamic allocation. */
6268 /* 33 = Pb in the data recuperation from block data */
6269 /* of coeff. of integration by GAUSS method. */
6270 /* >100 Pb in the evaluation of FONCNP, the returned error code */
6271 /* is equal to the error code of FONCNP + 100. */
6273 /* COMMONS USED : */
6274 /* ---------------- */
6276 /* REFERENCES CALLED : */
6277 /* ----------------------- */
6279 /* DESCRIPTION/NOTES/LIMITATIONS : */
6280 /* ----------------------------------- */
6281 /* --> The approximation part is done in the space of dimension */
6282 /* NDIMEN (the sum of dimensions of sub-spaces). For example : */
6283 /* If NBSESP=2 and NDIMSE(1)=3, NDIMSE(2)=2, there is smoothing with */
6284 /* NDIMEN=5. The result (in COURBE(NDIMEN,NCOEFF,i) ), will be */
6285 /* composed of the result of smoothing of 3D function in */
6286 /* COURBE(1:3,1:NCOEFF,i) and of smoothing of 2D function in */
6287 /* COURBE(4:5,1:NCOEFF,i). */
6289 /* --> Routine FONCNP should be declared EXTERNAL in the program */
6290 /* calling MMA2FNC. */
6292 /* --> Function FONCNP, declared externally, should be declared */
6293 /* IMPERATIVELY in form : */
6294 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIUOUV,TCONST,NBPTAB */
6295 /* ,TTABLE,IDERIU,IDERIV,IERCOD) */
6296 /* where the input arguments are : */
6297 /* - NDIMEN is integer defined as the sum of dimensions of */
6298 /* sub-spaces (i.e. total dimension of the problem). */
6299 /* - UINTFN(2) is a table of 2 reals containing the interval */
6300 /* by u where the function to be approximated is defined */
6301 /* (so it is equal to UIFONC). */
6302 /* - VINTFN(2) is a table of 2 reals containing the interval */
6303 /* by v where the function to be approximated is defined */
6304 /* (so it is equal to VIFONC). */
6305 /* - IIUOUV, shows that the points to be calculated have a constant U */
6306 /* (IIUOUV=1) or a constant V (IIUOUV=2). */
6307 /* - TCONST, real, value of the fixed discretisation parameter. Takes values */
6308 /* in (UINTFN(1),UINTFN(2)) if IIUOUV=1, */
6309 /* or in (VINTFN(1),VINTFN(2)) if IIUOUV=2. */
6310 /* - NBPTAB, the nb of point of discretisation following the free variable */
6311 /* : V if IIUOUV=1 or U if IIUOUV = 2. */
6312 /* - TTABLE, Table of NBPTAB parametres of discretisation. . */
6313 /* - IDERIU, integer, takes values between 0 (position) */
6314 /* and IORDREU (partial derivative of the function by u */
6315 /* of order IORDREU if IORDREU > 0). */
6316 /* - IDERIV, integer, takes values between 0 (position) */
6317 /* and IORDREV (partial derivative of the function by v */
6318 /* of order IORDREV if IORDREV > 0). */
6319 /* and the output arguments are : */
6320 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
6321 /* NBPTAB points calculated in FONCNP. */
6322 /* - IERCOD is, at output the error code of FONCNP. This code */
6323 /* (integer) should be strictly positive if there is a problem. */
6325 /* The input arguments SHOULD NOT BE modified under FONCNP.
6328 /* --> If IERCOD=-1, the required precision can't be reached (ERRMAX */
6329 /* is above EPSAPR on at least one sub-space), but
6331 /* one gives the best possible result for NCFLIM and EPSAPR */
6332 /* chosen by the user. In this case (and for IERCOD=0), there is a solution. */
6335 /* **********************************************************************
6337 /* Name of the routine */
6339 /* Parameter adjustments */
6344 errmoy_dim1 = *nbsesp;
6345 errmoy_offset = errmoy_dim1 + 1;
6346 errmoy -= errmoy_offset;
6347 errmax_dim1 = *nbsesp;
6348 errmax_offset = errmax_dim1 + 1;
6349 errmax -= errmax_offset;
6350 contr2_dim1 = *ndimen;
6351 contr2_dim2 = *iordre + 2;
6352 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
6353 contr2 -= contr2_offset;
6354 contr1_dim1 = *ndimen;
6355 contr1_dim2 = *iordre + 2;
6356 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
6357 contr1 -= contr1_offset;
6358 diftab_dim1 = *nbroot / 2 + 1;
6359 diftab_dim2 = *ndimen;
6360 diftab_offset = diftab_dim1 * (diftab_dim2 + 1);
6361 diftab -= diftab_offset;
6362 somtab_dim1 = *nbroot / 2 + 1;
6363 somtab_dim2 = *ndimen;
6364 somtab_offset = somtab_dim1 * (somtab_dim2 + 1);
6365 somtab -= somtab_offset;
6367 courbe_dim1 = *ncflim;
6368 courbe_dim2 = *ndimen;
6369 courbe_offset = courbe_dim1 * (courbe_dim2 + 1) + 1;
6370 courbe -= courbe_offset;
6371 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
6374 ibb = AdvApp2Var_SysBase::mnfndeb_();
6376 AdvApp2Var_SysBase::mgenmsg_("MMA2FNC", 7L);
6381 /* ---------------- Set to zero the coefficients of CURVE --------------
6384 ilong = *ndimen * *ncflim * *nbcrmx;
6385 AdvApp2Var_SysBase::mvriraz_(&ilong, &courbe[courbe_offset]);
6387 /* **********************************************************************
6389 /* -------------------------- Checking of entries ------------------
6391 /* **********************************************************************
6394 AdvApp2Var_MathBase::mmveps3_(&eps3);
6395 if ((d__1 = uvfonc[4] - uvfonc[3], advapp_abs(d__1)) < eps3) {
6398 if ((d__1 = uvfonc[6] - uvfonc[5], advapp_abs(d__1)) < eps3) {
6407 /* ********************************************************************** */
6408 /* ------------- Preparation of parameters of discretisation ----------- */
6409 /* **********************************************************************
6412 /* -- Allocation of a table of parameters and points of discretisation -- */
6413 /* --> For the parameters of discretisation. */
6415 /* --> For the points of discretisation in MMA1FDI and MMA1CDI and
6417 /* the auxiliary curve for MMAPCMP */
6418 ibid1 = *ndimen * (*nbroot + 2);
6419 ibid2 = ((*iordre + 1) << 1) * *nbroot;
6420 isz2 = advapp_max(ibid1,ibid2);
6421 ibid1 = (((*ncflim - 1) / 2 + 1) << 1) * *ndimen;
6422 isz2 = advapp_max(ibid1,isz2);
6423 /* --> To return the polynoms of hermit. */
6424 isz3 = ((*iordre + 1) << 2) * (*iordre + 1);
6425 /* --> For the Gauss coeff. of integration. */
6426 isz4 = (*nbroot / 2 + 1) * (*ndgjac + 1 - ((*iordre + 1) << 1));
6427 /* --> For the coeff of the curve in the base of Jacobi */
6428 isz5 = (*ndgjac + 1) * *ndimen;
6430 ndwrk = isz1 + isz2 + isz3 + isz4 + isz5;
6431 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6434 /* --> For the parameters of discretisation (NBROOT+2 extremities). */
6436 /* --> For the points of discretisation FPNTAB(NDIMEN,NBROOT+2), */
6437 /* FPNTAB(NBROOT,2*(IORDRE+1)) and for WRKAR of MMAPCMP. */
6439 /* --> For the polynoms of Hermit */
6441 /* --> For the Gauss coeff of integration. */
6443 /* --> For the curve in Jacobi. */
6446 /* ------------------ Initialisation of management of cuts ---------
6450 uvpav[0] = uvfonc[3];
6451 uvpav[1] = uvfonc[4];
6452 tabdec[0] = uvfonc[5];
6453 tabdec[1] = uvfonc[6];
6454 } else if (*isofav == 2) {
6455 tabdec[0] = uvfonc[3];
6456 tabdec[1] = uvfonc[4];
6457 uvpav[2] = uvfonc[5];
6458 uvpav[3] = uvfonc[6];
6466 /* **********************************************************************
6468 /* APPROXIMATION WITH CUTS */
6469 /* **********************************************************************
6473 /* --> When the top is reached, this is the end ! */
6474 if (nupil - *nbcrbe == 0) {
6479 uvpav[2] = tabdec[*nbcrbe];
6480 uvpav[3] = tabdec[*nbcrbe + 1];
6481 } else if (*isofav == 2) {
6482 uvpav[0] = tabdec[*nbcrbe];
6483 uvpav[1] = tabdec[*nbcrbe + 1];
6488 /* -------------------- Normalization of parameters -------------------- */
6490 mma1nop_(nbroot, &rootlg[1], uvpav, isofav, &wrkar[ipt1], &ier);
6495 /* -------------------- Discretisation of FONCNP ------------------------ */
6497 mma1fdi_(ndimen, uvpav, foncnp, isofav, tconst, nbroot, &wrkar[ipt1],
6498 iordre, ideriv, &wrkar[ipt2], &somtab[(ncb1 * somtab_dim2 + 1) *
6499 somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6500 contr1[(ncb1 * contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1
6501 * contr2_dim2 + 1) * contr2_dim1 + 1], iercod);
6506 /* -----------Cut the discretisation of constraints ------------*/
6509 mma1cdi_(ndimen, nbroot, &rootlg[1], iordre, &contr1[(ncb1 *
6510 contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1 *
6511 contr2_dim2 + 1) * contr2_dim1 + 1], &somtab[(ncb1 *
6512 somtab_dim2 + 1) * somtab_dim1], &diftab[(ncb1 * diftab_dim2
6513 + 1) * diftab_dim1], &wrkar[ipt2], &wrkar[ipt3], &ier);
6519 /* **********************************************************************
6521 /* -------------------- Calculate the curve of approximation -------------
6523 /* **********************************************************************
6526 mma1jak_(ndimen, nbroot, iordre, ndgjac, &somtab[(ncb1 * somtab_dim2 + 1)
6527 * somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6528 wrkar[ipt4], &wrkar[ipt5], &ier);
6533 /* **********************************************************************
6535 /* ---------------- Add polynom of interpolation -------------------
6537 /* **********************************************************************
6541 mma1cnt_(ndimen, iordre, &contr1[(ncb1 * contr1_dim2 + 1) *
6542 contr1_dim1 + 1], &contr2[(ncb1 * contr2_dim2 + 1) *
6543 contr2_dim1 + 1], &wrkar[ipt3], ndgjac, &wrkar[ipt5]);
6546 /* **********************************************************************
6548 /* --------------- Calculate Max and Average error ----------------------
6550 /* **********************************************************************
6553 mma1fer_(ndimen, nbsesp, &ndimse[1], iordre, ndgjac, &wrkar[ipt5], ncflim,
6554 &epsapr[1], &wrkar[ipt2], &errmax[ncb1 * errmax_dim1 + 1], &
6555 errmoy[ncb1 * errmoy_dim1 + 1], &ncoeff[ncb1], &ier);
6560 if (ier == 0 || (ier == -1 && nupil == *nbcrmx)) {
6562 /* ******************************************************************
6564 /* ----------------------- Compression du resultat ------------------
6566 /* ******************************************************************
6572 ncfja = *ndgjac + 1;
6573 /* -> Compression of result in WRKAR(IPT2) */
6576 AdvApp2Var_MathBase::mmapcmp_(ndimen,
6577 &ncfja, &ncoeff[ncb1], &wrkar[ipt5], &wrkar[ipt2]);
6579 AdvApp2Var_MathBase::mmapcmp_((integer*)ndimen,
6585 ilong = *ndimen * *ncflim;
6586 AdvApp2Var_SysBase::mvriraz_(&ilong, &wrkar[ipt5]);
6587 /* -> Passage to canonic base (-1,1) (result in WRKAR(IPT5)).
6589 ndgre = ncoeff[ncb1] - 1;
6591 for (nd = 1; nd <= i__1; ++nd) {
6592 iptt = ipt2 + ((nd - 1) << 1) * (ndgre / 2 + 1);
6593 jptt = ipt5 + (nd - 1) * ncoeff[ncb1];
6594 AdvApp2Var_MathBase::mmjacan_(iordre, &ndgre, &wrkar[iptt], &wrkar[jptt]);
6598 /* -> Store the calculated curve */
6600 AdvApp2Var_MathBase::mmfmca8_(&ncoeff[ncb1], ndimen, &ibid1, ncflim, ndimen, &ibid1, &
6601 wrkar[ipt5], &courbe[(ncb1 * courbe_dim2 + 1) * courbe_dim1 +
6604 /* -> Before normalization of constraints on (-1,1), recalculate */
6605 /* the true constraints. */
6607 for (ii = 0; ii <= i__1; ++ii) {
6608 mma1noc_(uv11, ndimen, &ii, &contr1[(ii + 1 + ncb1 * contr1_dim2)
6609 * contr1_dim1 + 1], uvpav, isofav, ideriv, &contr1[(ii +
6610 1 + ncb1 * contr1_dim2) * contr1_dim1 + 1]);
6611 mma1noc_(uv11, ndimen, &ii, &contr2[(ii + 1 + ncb1 * contr2_dim2)
6612 * contr2_dim1 + 1], uvpav, isofav, ideriv, &contr2[(ii +
6613 1 + ncb1 * contr2_dim2) * contr2_dim1 + 1]);
6617 ibid1 = (*nbroot / 2 + 1) * *ndimen;
6618 mma1noc_(uv11, &ibid1, &ii, &somtab[(ncb1 * somtab_dim2 + 1) *
6619 somtab_dim1], uvpav, isofav, ideriv, &somtab[(ncb1 *
6620 somtab_dim2 + 1) * somtab_dim1]);
6621 mma1noc_(uv11, &ibid1, &ii, &diftab[(ncb1 * diftab_dim2 + 1) *
6622 diftab_dim1], uvpav, isofav, ideriv, &diftab[(ncb1 *
6623 diftab_dim2 + 1) * diftab_dim1]);
6626 for (nd = 1; nd <= i__1; ++nd) {
6627 mma1noc_(uv11, &ncoeff[ncb1], &ii, &courbe[(nd + ncb1 *
6628 courbe_dim2) * courbe_dim1 + 1], uvpav, isofav, ideriv, &
6629 courbe[(nd + ncb1 * courbe_dim2) * courbe_dim1 + 1]);
6633 /* -> Update the nb of already created curves */
6636 /* -> ...otherwise try to cut the current interval in 2... */
6638 tmil = (tabdec[*nbcrbe + 1] + tabdec[*nbcrbe]) / 2.;
6641 ilong = (nupil - *nbcrbe) << 3;
6642 AdvApp2Var_SysBase::mcrfill_(&ilong, &tabdec[ideb],&tabdec[ideb1]);
6643 tabdec[ideb] = tmil;
6647 /* ---------- Make approximation of the rest -----------
6652 /* --------------------- Return code of error -----------------
6654 /* --> Pb with dynamic allocation */
6658 /* --> Inputs incoherent. */
6663 /* -------------------------- Dynamic desallocation -------------------
6668 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6675 /* ------------------------------ The end -------------------------------
6680 AdvApp2Var_SysBase::maermsg_("MMA2FNC", iercod, 7L);
6683 AdvApp2Var_SysBase::mgsomsg_("MMA2FNC", 7L);
6688 //=======================================================================
6689 //function : mma2fx6_
6691 //=======================================================================
6692 int AdvApp2Var_ApproxF2var::mma2fx6_(integer *ncfmxu,
6709 /* System generated locals */
6710 integer epsfro_dim1, epsfro_offset, patcan_dim1, patcan_dim2, patcan_dim3,
6711 patcan_dim4, patcan_offset, errmax_dim1, errmax_dim2,
6712 errmax_offset, ncoefu_dim1, ncoefu_offset, ncoefv_dim1,
6713 ncoefv_offset, i__1, i__2, i__3, i__4, i__5;
6714 doublereal d__1, d__2;
6716 /* Local variables */
6717 integer idim, ncfu, ncfv, id, ii, nd, jj, ku, kv, ns, ibb;
6721 /* **********************************************************************
6726 /* Reduction of degree when the squares are the squares of constraints. */
6730 /* TOUS,AB_SPECIFI::CARREAU&,REDUCTION,&CARREAU */
6732 /* INPUT ARGUMENTS : */
6733 /* ------------------ */
6734 /* NCFMXU: Max Nb of coeff by u of solution P(u,v) (table */
6735 /* PATCAN). This argument serves only to declare the size of this table. */
6736 /* NCFMXV: Max Nb of coeff by v of solution P(u,v) (table */
6737 /* PATCAN). This argument serves only to declare the size of this table. */
6738 /* NDIMEN: Total dimension of the space where the processed function */
6739 /* takes its values.(sum of dimensions of sub-spaces) */
6740 /* NBSESP: Nb of independent sub-spaces where the errors are measured. */
6741 /* NDIMSE: Table of dimensions of NBSESP sub-spaces. */
6742 /* NBUPAT: Nb of square solution by u. */
6743 /* NBVPAT: Nb of square solution by v. */
6744 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
6745 /* = 0, the extremities of iso-V are calculated */
6746 /* = 1, additionally the 1st derivative in the direction of iso-V is calculated */
6747 /* = 2, additionally the 2nd derivative in the direction of iso-V is calculated */
6748 /* IORDRV: Ordre de contrainte impose aux extremites de l'iso-U */
6749 /* = 0, on calcule les extremites de l'iso-U. */
6750 /* = 1, additionally the 1st derivative in the direction of iso-U is calculated */
6751 /* = 2, additionally the 2nd derivative in the direction of iso-U is calculated */
6752 /* EPSAPR: Table of imposed precisions, sub-space by sub-space. */
6753 /* EPSFRO: Table of imposed precisions, sub-space by sub-space on the limits of squares. */
6754 /* PATCAN: Table of coeff. in the canonic base of squares P(u,v) calculated for (u,v) in (-1,1). */
6755 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6756 /* committed in the approximation of F(u,v) by P(u,v). */
6757 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6758 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6760 /* OUTPUT ARGUMENTS : */
6761 /* ------------------- */
6762 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6763 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6765 /* COMMONS USED : */
6766 /* ---------------- */
6768 /* REFERENCES CALLED : */
6769 /* --------------------- */
6771 /* DESCRIPTION/NOTES/LIMITATIONS : */
6772 /* ------------------------------- */
6774 /* **********************************************************************
6777 /* Name of the routine */
6780 /* Parameter adjustments */
6781 epsfro_dim1 = *nbsesp;
6782 epsfro_offset = epsfro_dim1 * 5 + 1;
6783 epsfro -= epsfro_offset;
6786 ncoefv_dim1 = *nbupat;
6787 ncoefv_offset = ncoefv_dim1 + 1;
6788 ncoefv -= ncoefv_offset;
6789 ncoefu_dim1 = *nbupat;
6790 ncoefu_offset = ncoefu_dim1 + 1;
6791 ncoefu -= ncoefu_offset;
6792 errmax_dim1 = *nbsesp;
6793 errmax_dim2 = *nbupat;
6794 errmax_offset = errmax_dim1 * (errmax_dim2 + 1) + 1;
6795 errmax -= errmax_offset;
6796 patcan_dim1 = *ncfmxu;
6797 patcan_dim2 = *ncfmxv;
6798 patcan_dim3 = *ndimen;
6799 patcan_dim4 = *nbupat;
6800 patcan_offset = patcan_dim1 * (patcan_dim2 * (patcan_dim3 * (patcan_dim4
6802 patcan -= patcan_offset;
6805 ibb = AdvApp2Var_SysBase::mnfndeb_();
6807 AdvApp2Var_SysBase::mgenmsg_("MMA2FX6", 7L);
6812 for (jj = 1; jj <= i__1; ++jj) {
6814 for (ii = 1; ii <= i__2; ++ii) {
6815 ncfu = ncoefu[ii + jj * ncoefu_dim1];
6816 ncfv = ncoefv[ii + jj * ncoefv_dim1];
6818 /* ********************************************************************** */
6819 /* -------------------- Reduction of degree by U ------------------------- */
6820 /* ********************************************************************** */
6823 if (ncfu <= (*iordru + 1) << 1 && ncfu > 2) {
6827 for (ns = 1; ns <= i__3; ++ns) {
6830 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6831 tol = advapp_min(d__1,d__2);
6833 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6834 tol = advapp_min(d__1,d__2);
6836 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6837 tol = advapp_min(d__1,d__2);
6839 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6840 tol = advapp_min(d__1,d__2);
6841 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6844 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6845 tol = advapp_min(d__1,d__2);
6847 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6848 tol = advapp_min(d__1,d__2);
6850 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6851 tol = advapp_min(d__1,d__2);
6853 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6854 tol = advapp_min(d__1,d__2);
6859 for (nd = 1; nd <= i__4; ++nd) {
6862 for (kv = 1; kv <= i__5; ++kv) {
6863 bid += (d__1 = patcan[ncfu + (kv + (id + (ii + jj
6864 * patcan_dim4) * patcan_dim3) *
6865 patcan_dim2) * patcan_dim1], advapp_abs(d__1));
6871 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6872 errmax_dim2) * errmax_dim1]) {
6883 /* ********************************************************************** */
6884 /* -------------------- Reduction of degree by V ------------------------- */
6885 /* ********************************************************************** */
6888 if (ncfv <= (*iordrv + 1) << 1 && ncfv > 2) {
6892 for (ns = 1; ns <= i__3; ++ns) {
6895 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6896 tol = advapp_min(d__1,d__2);
6898 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6899 tol = advapp_min(d__1,d__2);
6901 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6902 tol = advapp_min(d__1,d__2);
6904 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6905 tol = advapp_min(d__1,d__2);
6906 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6909 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6910 tol = advapp_min(d__1,d__2);
6912 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6913 tol = advapp_min(d__1,d__2);
6915 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6916 tol = advapp_min(d__1,d__2);
6918 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6919 tol = advapp_min(d__1,d__2);
6924 for (nd = 1; nd <= i__4; ++nd) {
6927 for (ku = 1; ku <= i__5; ++ku) {
6928 bid += (d__1 = patcan[ku + (ncfv + (id + (ii + jj
6929 * patcan_dim4) * patcan_dim3) *
6930 patcan_dim2) * patcan_dim1], advapp_abs(d__1));
6936 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6937 errmax_dim2) * errmax_dim1]) {
6948 /* --- Return the nbs of coeff. and pass to the next square --- */
6951 ncoefu[ii + jj * ncoefu_dim1] = advapp_max(ncfu,2);
6952 ncoefv[ii + jj * ncoefv_dim1] = advapp_max(ncfv,2);
6958 /* ------------------------------ The End -------------------------------
6962 AdvApp2Var_SysBase::mgsomsg_("MMA2FX6", 7L);
6968 //=======================================================================
6969 //function : mma2jmx_
6971 //=======================================================================
6972 int AdvApp2Var_ApproxF2var::mma2jmx_(integer *ndgjac,
6976 /* Initialized data */
6978 static doublereal xmax2[57] = { .9682458365518542212948163499456,
6979 .986013297183269340427888048593603,
6980 1.07810420343739860362585159028115,
6981 1.17325804490920057010925920756025,
6982 1.26476561266905634732910520370741,
6983 1.35169950227289626684434056681946,
6984 1.43424378958284137759129885012494,
6985 1.51281316274895465689402798226634,
6986 1.5878364329591908800533936587012,
6987 1.65970112228228167018443636171226,
6988 1.72874345388622461848433443013543,
6989 1.7952515611463877544077632304216,
6990 1.85947199025328260370244491818047,
6991 1.92161634324190018916351663207101,
6992 1.98186713586472025397859895825157,
6993 2.04038269834980146276967984252188,
6994 2.09730119173852573441223706382076,
6995 2.15274387655763462685970799663412,
6996 2.20681777186342079455059961912859,
6997 2.25961782459354604684402726624239,
6998 2.31122868752403808176824020121524,
6999 2.36172618435386566570998793688131,
7000 2.41117852396114589446497298177554,
7001 2.45964731268663657873849811095449,
7002 2.50718840313973523778244737914028,
7003 2.55385260994795361951813645784034,
7004 2.59968631659221867834697883938297,
7005 2.64473199258285846332860663371298,
7006 2.68902863641518586789566216064557,
7007 2.73261215675199397407027673053895,
7008 2.77551570192374483822124304745691,
7009 2.8177699459714315371037628127545,
7010 2.85940333797200948896046563785957,
7011 2.90044232019793636101516293333324,
7012 2.94091151970640874812265419871976,
7013 2.98083391718088702956696303389061,
7014 3.02023099621926980436221568258656,
7015 3.05912287574998661724731962377847,
7016 3.09752842783622025614245706196447,
7017 3.13546538278134559341444834866301,
7018 3.17295042316122606504398054547289,
7019 3.2099992681699613513775259670214,
7020 3.24662674946606137764916854570219,
7021 3.28284687953866689817670991319787,
7022 3.31867291347259485044591136879087,
7023 3.35411740487202127264475726990106,
7024 3.38919225660177218727305224515862,
7025 3.42390876691942143189170489271753,
7026 3.45827767149820230182596660024454,
7027 3.49230918177808483937957161007792,
7028 3.5260130200285724149540352829756,
7029 3.55939845146044235497103883695448,
7030 3.59247431368364585025958062194665,
7031 3.62524904377393592090180712976368,
7032 3.65773070318071087226169680450936,
7033 3.68992700068237648299565823810245,
7034 3.72184531357268220291630708234186 };
7035 static doublereal xmax4[55] = { 1.1092649593311780079813740546678,
7036 1.05299572648705464724876659688996,
7037 1.0949715351434178709281698645813,
7038 1.15078388379719068145021100764647,
7039 1.2094863084718701596278219811869,
7040 1.26806623151369531323304177532868,
7041 1.32549784426476978866302826176202,
7042 1.38142537365039019558329304432581,
7043 1.43575531950773585146867625840552,
7044 1.48850442653629641402403231015299,
7045 1.53973611681876234549146350844736,
7046 1.58953193485272191557448229046492,
7047 1.63797820416306624705258190017418,
7048 1.68515974143594899185621942934906,
7049 1.73115699602477936547107755854868,
7050 1.77604489805513552087086912113251,
7051 1.81989256661534438347398400420601,
7052 1.86276344480103110090865609776681,
7053 1.90471563564740808542244678597105,
7054 1.94580231994751044968731427898046,
7055 1.98607219357764450634552790950067,
7056 2.02556989246317857340333585562678,
7057 2.06433638992049685189059517340452,
7058 2.10240936014742726236706004607473,
7059 2.13982350649113222745523925190532,
7060 2.17661085564771614285379929798896,
7061 2.21280102016879766322589373557048,
7062 2.2484214321456956597803794333791,
7063 2.28349755104077956674135810027654,
7064 2.31805304852593774867640120860446,
7065 2.35210997297725685169643559615022,
7066 2.38568889602346315560143377261814,
7067 2.41880904328694215730192284109322,
7068 2.45148841120796359750021227795539,
7069 2.48374387161372199992570528025315,
7070 2.5155912654873773953959098501893,
7071 2.54704548720896557684101746505398,
7072 2.57812056037881628390134077704127,
7073 2.60882970619319538196517982945269,
7074 2.63918540521920497868347679257107,
7075 2.66919945330942891495458446613851,
7076 2.69888301230439621709803756505788,
7077 2.72824665609081486737132853370048,
7078 2.75730041251405791603760003778285,
7079 2.78605380158311346185098508516203,
7080 2.81451587035387403267676338931454,
7081 2.84269522483114290814009184272637,
7082 2.87060005919012917988363332454033,
7083 2.89823818258367657739520912946934,
7084 2.92561704377132528239806135133273,
7085 2.95274375377994262301217318010209,
7086 2.97962510678256471794289060402033,
7087 3.00626759936182712291041810228171,
7088 3.03267744830655121818899164295959,
7089 3.05886060707437081434964933864149 };
7090 static doublereal xmax6[53] = { 1.21091229812484768570102219548814,
7091 1.11626917091567929907256116528817,
7092 1.1327140810290884106278510474203,
7093 1.1679452722668028753522098022171,
7094 1.20910611986279066645602153641334,
7095 1.25228283758701572089625983127043,
7096 1.29591971597287895911380446311508,
7097 1.3393138157481884258308028584917,
7098 1.3821288728999671920677617491385,
7099 1.42420414683357356104823573391816,
7100 1.46546895108549501306970087318319,
7101 1.50590085198398789708599726315869,
7102 1.54550385142820987194251585145013,
7103 1.58429644271680300005206185490937,
7104 1.62230484071440103826322971668038,
7105 1.65955905239130512405565733793667,
7106 1.69609056468292429853775667485212,
7107 1.73193098017228915881592458573809,
7108 1.7671112206990325429863426635397,
7109 1.80166107681586964987277458875667,
7110 1.83560897003644959204940535551721,
7111 1.86898184653271388435058371983316,
7112 1.90180515174518670797686768515502,
7113 1.93410285411785808749237200054739,
7114 1.96589749778987993293150856865539,
7115 1.99721027139062501070081653790635,
7116 2.02806108474738744005306947877164,
7117 2.05846864831762572089033752595401,
7118 2.08845055210580131460156962214748,
7119 2.11802334209486194329576724042253,
7120 2.14720259305166593214642386780469,
7121 2.17600297710595096918495785742803,
7122 2.20443832785205516555772788192013,
7123 2.2325216999457379530416998244706,
7124 2.2602654243075083168599953074345,
7125 2.28768115912702794202525264301585,
7126 2.3147799369092684021274946755348,
7127 2.34157220782483457076721300512406,
7128 2.36806787963276257263034969490066,
7129 2.39427635443992520016789041085844,
7130 2.42020656255081863955040620243062,
7131 2.44586699364757383088888037359254,
7132 2.47126572552427660024678584642791,
7133 2.49641045058324178349347438430311,
7134 2.52130850028451113942299097584818,
7135 2.54596686772399937214920135190177,
7136 2.5703922285006754089328998222275,
7137 2.59459096001908861492582631591134,
7138 2.61856915936049852435394597597773,
7139 2.64233265984385295286445444361827,
7140 2.66588704638685848486056711408168,
7141 2.68923766976735295746679957665724,
7142 2.71238965987606292679677228666411 };
7144 /* System generated locals */
7147 /* Local variables */
7153 /* **********************************************************************
7158 /* Calculate the max of Jacobo polynoms multiplied by the weight on */
7159 /* (-1,1) for order 0,4,6 or Legendre. */
7163 /* LEGENDRE,APPROXIMATION,ERREUR. */
7165 /* INPUT ARGUMENTS : */
7166 /* ------------------ */
7167 /* NDGJAC: Nb of Jacobi coeff. of approximation. */
7168 /* IORDRE: Order of continuity (from -1 to 2) */
7170 /* OUTPUT ARGUMENTS : */
7171 /* ------------------- */
7172 /* XJACMX: Table of maximums of Jacobi polynoms. */
7174 /* COMMONS USED : */
7175 /* ---------------- */
7177 /* REFERENCES CALLED : */
7178 /* --------------------- */
7180 /* DESCRIPTION/NOTES/LIMITATIONS : */
7181 /* ----------------------------------- */
7184 /* ***********************************************************************
7186 /* Name of the routine */
7187 /* ----------------------------- Initialisations ------------------------
7190 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7192 AdvApp2Var_SysBase::mgenmsg_("MMA2JMX", 7L);
7195 numax = *ndgjac - ((*iordre + 1) << 1);
7196 if (*iordre == -1) {
7198 for (ii = 0; ii <= i__1; ++ii) {
7199 bid = (ii * 2. + 1.) / 2.;
7200 xjacmx[ii] = sqrt(bid);
7203 } else if (*iordre == 0) {
7205 for (ii = 0; ii <= i__1; ++ii) {
7206 xjacmx[ii] = xmax2[ii];
7209 } else if (*iordre == 1) {
7211 for (ii = 0; ii <= i__1; ++ii) {
7212 xjacmx[ii] = xmax4[ii];
7215 } else if (*iordre == 2) {
7217 for (ii = 0; ii <= i__1; ++ii) {
7218 xjacmx[ii] = xmax6[ii];
7223 /* ------------------------- The end ------------------------------------
7227 AdvApp2Var_SysBase::mgsomsg_("MMA2JMX", 7L);
7232 //=======================================================================
7233 //function : mma2moy_
7235 //=======================================================================
7236 int mma2moy_(integer *ndgumx,
7248 /* System generated locals */
7249 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
7251 /* Local variables */
7253 integer minu, minv, idebu, idebv, ii, nd, jj;
7254 doublereal bid0, bid1;
7257 /* **********************************************************************
7262 /* Calculate the average approximation error made when only */
7263 /* the coefficients of PATJAC of degree between */
7264 /* 2*(IORDRU+1) and MINDGU by U and 2*(IORDRV+1) and MINDGV by V are preserved. */
7268 /* LEGENDRE,APPROXIMATION, AVERAGE ERROR */
7270 /* INPUT ARGUMENTS : */
7271 /* ------------------ */
7272 /* NDGUMX: Dimension by U of table PATJAC. */
7273 /* NDGVMX: Dimension by V of table PATJAC. */
7274 /* NDIMEN: Dimension of the space. */
7275 /* MINDGU: Lower limit of the index by U of PATJAC coeff to be taken into account. */
7276 /* MAXDGU: Upper limit of the index by U of PATJAC coeff to be taken into account. */
7277 /* MINDGV: Lower limit of the index by V of PATJAC coeff to be taken into account. */
7278 /* MAXDGV: Upper limit of the index by V of PATJAC coeff to be taken into account. */
7279 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
7280 /* IORDRV: Order of continuity by V provided by square PATJAC (from -1 to 2) */
7281 /* PATJAC: Table of coeff. of the approximation square with */
7282 /* constraints of order IORDRU by U and IORDRV by V. */
7284 /* OUTPUT ARGUMENTS : */
7285 /* ------------------- */
7286 /* ERRMOY: Average error commited by preserving only the coeff of */
7287 /* PATJAC 2*(IORDRU+1) in MINDGU by U and 2*(IORDRV+1) in MINDGV by V. */
7289 /* COMMONS USED : */
7290 /* ---------------- */
7292 /* REFERENCES CALLED : */
7293 /* --------------------- */
7295 /* DESCRIPTION/NOTES/LIMITATIONS : */
7296 /* ----------------------------------- */
7297 /* Table PATJAC stores the coeff. Cij of */
7298 /* approximation square F(U,V). Indexes i and j show the degree by */
7299 /* U and by V of the base polynoms. These base polynoms are in the form: */
7301 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
7303 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
7304 /* IORDRU+1 (the same by V by replacing U by V in the above expression). */
7306 /* The contribution to the average error of term Cij when */
7307 /* it is removed from PATJAC is Cij*Cij. */
7310 /* ***********************************************************************
7312 /* Name of the routine */
7315 /* ----------------------------- Initialisations ------------------------
7318 /* Parameter adjustments */
7319 patjac_dim1 = *ndgumx + 1;
7320 patjac_dim2 = *ndgvmx + 1;
7321 patjac_offset = patjac_dim1 * patjac_dim2;
7322 patjac -= patjac_offset;
7325 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7327 AdvApp2Var_SysBase::mgenmsg_("MMA2MOY", 7L);
7330 idebu = (*iordru + 1) << 1;
7331 idebv = (*iordrv + 1) << 1;
7332 minu = advapp_max(idebu,*mindgu);
7333 minv = advapp_max(idebv,*mindgv);
7337 /* ------------------ Calculation of the upper bound of the average error ------------ */
7338 /* -------------------- when the coeff. of indexes from MINDGU to MAXDGU ------ */
7339 /* ---------------- by U and of indexes from MINDGV to MAXDGV by V are removed -------------- */
7342 for (nd = 1; nd <= i__1; ++nd) {
7344 for (jj = minv; jj <= i__2; ++jj) {
7346 for (ii = idebu; ii <= i__3; ++ii) {
7347 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7348 bid0 += bid1 * bid1;
7357 for (nd = 1; nd <= i__1; ++nd) {
7359 for (jj = idebv; jj <= i__2; ++jj) {
7361 for (ii = minu; ii <= i__3; ++ii) {
7362 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7363 bid0 += bid1 * bid1;
7371 /* ----------------------- Calculation of the average error -------------
7375 *errmoy = sqrt(bid0);
7377 /* ------------------------- The end ------------------------------------
7381 AdvApp2Var_SysBase::mgsomsg_("MMA2MOY", 7L);
7386 //=======================================================================
7387 //function : mma2roo_
7389 //=======================================================================
7390 int AdvApp2Var_ApproxF2var::mma2roo_(integer *nbpntu,
7395 /* System generated locals */
7398 /* Local variables */
7401 /* **********************************************************************
7406 /* Return roots of Legendre for discretisations. */
7410 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
7412 /* INPUT ARGUMENTS : */
7413 /* ------------------ */
7414 /* NBPNTU: Nb of INTERNAL parameters of discretization BY U. */
7415 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7416 /* NBPNTV: Nb of INTERNAL parameters of discretization BY V. */
7417 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7419 /* OUTPUT ARGUMENTS : */
7420 /* ------------------- */
7421 /* UROOTL: Table of parameters of discretisation ON (-1,1) BY U.
7423 /* VROOTL: Table of parameters of discretisation ON (-1,1) BY V.
7426 /* COMMONS USED : */
7427 /* ---------------- */
7429 /* REFERENCES CALLED : */
7430 /* --------------------- */
7432 /* DESCRIPTION/NOTES/LIMITATIONS : */
7433 /* ----------------------------------- */
7436 /* **********************************************************************
7439 /* Name of the routine */
7442 /* Parameter adjustments */
7447 ibb = AdvApp2Var_SysBase::mnfndeb_();
7449 AdvApp2Var_SysBase::mgenmsg_("MMA2ROO", 7L);
7452 /* ---------------- Return the POSITIVE roots on U ------------------
7455 AdvApp2Var_MathBase::mmrtptt_(nbpntu, &urootl[(*nbpntu + 1) / 2 + 1]);
7457 for (ii = 1; ii <= i__1; ++ii) {
7458 urootl[ii] = -urootl[*nbpntu - ii + 1];
7461 if (*nbpntu % 2 == 1) {
7462 urootl[*nbpntu / 2 + 1] = 0.;
7465 /* ---------------- Return the POSITIVE roots on V ------------------
7468 AdvApp2Var_MathBase::mmrtptt_(nbpntv, &vrootl[(*nbpntv + 1) / 2 + 1]);
7470 for (ii = 1; ii <= i__1; ++ii) {
7471 vrootl[ii] = -vrootl[*nbpntv - ii + 1];
7474 if (*nbpntv % 2 == 1) {
7475 vrootl[*nbpntv / 2 + 1] = 0.;
7478 /* ------------------------------ The End -------------------------------
7482 AdvApp2Var_SysBase::mgsomsg_("MMA2ROO", 7L);
7486 //=======================================================================
7487 //function : mmmapcoe_
7489 //=======================================================================
7490 int mmmapcoe_(integer *ndim,
7500 /* System generated locals */
7501 integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
7502 crvjac_dim1, crvjac_offset, gsstab_dim1, i__1, i__2, i__3;
7504 /* Local variables */
7505 integer igss, ikdeb;
7507 integer nd, ik, ir, nbroot, ibb;
7509 /* **********************************************************************
7514 /* Calculate the coefficients of polinomial approximation curve */
7515 /* of degree NDGJAC by the method of smallest squares starting from */
7516 /* the discretization of function on the roots of Legendre polynom */
7517 /* of degree NBPNTS. */
7521 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
7523 /* INPUT ARGUMENTS : */
7524 /* ------------------ */
7525 /* NDIM : Dimension of the space. */
7526 /* NDGJAC : Max Degree of the polynom of approximation. */
7527 /* The representation in the orthogonal base starts from degree */
7528 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7529 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7530 /* IORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7531 /* to step of constraints, C0,C1 or C2. */
7532 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7533 /* are calculated the coefficients of integration by */
7534 /* Gauss method. It is required to set NBPNTS=30,40,50 or 61 */
7535 /* and NDGJAC < NBPNTS. */
7536 /* SOMTAB : Table of F(ti)+F(-ti) with ti in ROOTAB. */
7537 /* DIFTAB : Table of F(ti)-F(-ti) with ti in ROOTAB. */
7538 /* GSSTAB(i,k) : Table of coefficients of integration by the Gauss method : */
7539 /* i varies from 0 to NBPNTS and */
7540 /* k varies from 0 to NDGJAC-2*(JORDRE+1). */
7542 /* OUTPUT ARGUMENTSE : */
7543 /* ------------------- */
7544 /* CRVJAC : Curve of approximation of FONCNP with eventually */
7545 /* taking into account of constraints at the extremities. */
7546 /* This curve is of degree NDGJAC. */
7548 /* COMMONS USED : */
7549 /* ---------------- */
7551 /* REFERENCES CALLED : */
7552 /* --------------------- */
7554 /* DESCRIPTION/NOTES/LIMITATIONS : */
7555 /* ------------------------------- */
7557 /* **********************************************************************
7560 /* Name of the routine */
7562 /* Parameter adjustments */
7563 crvjac_dim1 = *ndgjac + 1;
7564 crvjac_offset = crvjac_dim1;
7565 crvjac -= crvjac_offset;
7566 gsstab_dim1 = *nbpnts / 2 + 1;
7567 diftab_dim1 = *nbpnts / 2 + 1;
7568 diftab_offset = diftab_dim1;
7569 diftab -= diftab_offset;
7570 somtab_dim1 = *nbpnts / 2 + 1;
7571 somtab_offset = somtab_dim1;
7572 somtab -= somtab_offset;
7575 ibb = AdvApp2Var_SysBase::mnfndeb_();
7577 AdvApp2Var_SysBase::mgenmsg_("MMMAPCO", 7L);
7579 ikdeb = (*iordre + 1) << 1;
7580 nbroot = *nbpnts / 2;
7583 for (nd = 1; nd <= i__1; ++nd) {
7585 /* ----------------- Calculate the coefficients of even degree ----------
7589 for (ik = ikdeb; ik <= i__2; ik += 2) {
7593 for (ir = 1; ir <= i__3; ++ir) {
7594 bidon += somtab[ir + nd * somtab_dim1] * gsstab[ir + igss *
7598 crvjac[ik + nd * crvjac_dim1] = bidon;
7602 /* --------------- Calculate the coefficients of uneven degree ----------
7606 for (ik = ikdeb + 1; ik <= i__2; ik += 2) {
7610 for (ir = 1; ir <= i__3; ++ir) {
7611 bidon += diftab[ir + nd * diftab_dim1] * gsstab[ir + igss *
7615 crvjac[ik + nd * crvjac_dim1] = bidon;
7622 /* ------- Add terms connected to the supplementary root (0.D0) ------ */
7623 /* ----------- of Legendre polynom of uneven degree NBPNTS -----------
7626 if (*nbpnts % 2 == 0) {
7630 for (nd = 1; nd <= i__1; ++nd) {
7632 for (ik = ikdeb; ik <= i__2; ik += 2) {
7634 crvjac[ik + nd * crvjac_dim1] += somtab[nd * somtab_dim1] *
7635 gsstab[igss * gsstab_dim1];
7641 /* ------------------------------ The end -------------------------------
7646 AdvApp2Var_SysBase::mgsomsg_("MMMAPCO", 7L);
7650 //=======================================================================
7651 //function : mmaperm_
7653 //=======================================================================
7654 int mmaperm_(integer *ncofmx,
7662 /* System generated locals */
7663 integer crvjac_dim1, crvjac_offset, i__1, i__2;
7665 /* Local variables */
7667 integer i__, ia, nd, ncfcut, ibb;
7670 /* **********************************************************************
7675 /* Calculate the square root of the average quadratic error */
7676 /* of approximation done when only the */
7677 /* first NCFNEW coefficients of a curve of degree NCOEFF-1 */
7678 /* written in NORMALIZED Jacobi base of order 2*(IORDRE+1) are preserved. */
7682 /* LEGENDRE,POLYGONE,APPROXIMATION,ERREUR. */
7684 /* INPUT ARGUMENTS : */
7685 /* ------------------ */
7686 /* NCOFMX : Maximum degree of the curve. */
7687 /* NDIM : Dimension of the space. */
7688 /* NCOEFF : Degree +1 of the curve. */
7689 /* IORDRE : Order of constraint of continuity at the extremities. */
7690 /* CRVJAC : The curve the degree which of will be lowered. */
7691 /* NCFNEW : Degree +1 of the resulting polynom. */
7693 /* OUTPUT ARGUMENTS : */
7694 /* ------------------- */
7695 /* ERRMOY : Average precision of approximation. */
7697 /* COMMONS USED : */
7698 /* ---------------- */
7700 /* REFERENCES CALLED : */
7701 /* ----------------------- */
7703 /* DESCRIPTION/NOTES/LIMITATIONS : */
7704 /* ----------------------------------- */
7706 /* ***********************************************************************
7709 /* Name of the routine */
7711 /* Parameter adjustments */
7712 crvjac_dim1 = *ncofmx;
7713 crvjac_offset = crvjac_dim1 + 1;
7714 crvjac -= crvjac_offset;
7717 ibb = AdvApp2Var_SysBase::mnfndeb_();
7719 AdvApp2Var_SysBase::mgenmsg_("MMAPERM", 7L);
7722 /* --------- Minimum degree that can be reached : Stop at 1 or IA -------
7725 ia = (*iordre + 1) << 1;
7727 if (*ncfnew + 1 > ncfcut) {
7728 ncfcut = *ncfnew + 1;
7731 /* -------------- Elimination of coefficients of high degree ------------ */
7732 /* ----------- Loop on the series of Jacobi :NCFCUT --> NCOEFF --------- */
7737 for (nd = 1; nd <= i__1; ++nd) {
7739 for (i__ = ncfcut; i__ <= i__2; ++i__) {
7740 bidj = crvjac[i__ + nd * crvjac_dim1];
7747 /* ----------- Square Root of average quadratic error e -----------
7751 *errmoy = sqrt(bid);
7753 /* ------------------------------- The end ------------------------------
7757 AdvApp2Var_SysBase::mgsomsg_("MMAPERM", 7L);
7761 //=======================================================================
7762 //function : mmapptt_
7764 //=======================================================================
7765 int AdvApp2Var_ApproxF2var::mmapptt_(const integer *ndgjac,
7766 const integer *nbpnts,
7767 const integer *jordre,
7771 /* System generated locals */
7772 integer cgauss_dim1, i__1;
7774 /* Local variables */
7775 integer kjac, iptt, ipdb0, infdg, iptdb, mxjac, ilong, ibb;
7777 /* **********************************************************************
7782 /* Load the elements required for integration by */
7783 /* Gauss method to obtain the coefficients in the base of */
7784 /* Legendre of the approximation by the least squares of a */
7785 /* function. The elements are stored in commons MMAPGSS */
7786 /* (case without constraint), MMAPGS0 (constraints C0), MMAPGS1 */
7787 /* (constraints C1) and MMAPGS2 (constraints C2). */
7791 /* INTEGRATION,GAUSS,JACOBI */
7793 /* INPUT ARGUMENTS : */
7794 /* ------------------ */
7795 /* NDGJAC : Max degree of the polynom of approximation. */
7796 /* The representation in orthogonal base goes from degree */
7797 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7798 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7799 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7800 /* are calculated the coefficients of integration by the */
7801 /* method of Gauss. It is required that NBPNTS=8,10,15,20,25, */
7802 /* 30,40,50 or 61 and NDGJAC < NBPNTS. */
7803 /* JORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7804 /* to step of constraints C0,C1 or C2. */
7806 /* OUTPUT ARGUMENTS : */
7807 /* ------------------- */
7808 /* CGAUSS(i,k) : Table of coefficients of integration by */
7809 /* Gauss method : i varies from 0 to the integer part */
7810 /* of NBPNTS/2 and k varies from 0 to NDGJAC-2*(JORDRE+1). */
7811 /* These are the coeff. of integration associated to */
7812 /* positive roots of the polynom of Legendre of degree */
7813 /* NBPNTS. CGAUSS(0,k) contains coeff. */
7814 /* of integration associated to root t = 0 when */
7815 /* NBPNTS is uneven. */
7816 /* IERCOD : Error code. */
7818 /* = 11 NBPNTS is not 8,10,15,20,25,30,40,50 or 61. */
7819 /* = 21 JORDRE is not -1,0,1 or 2. */
7820 /* = 31 NDGJAC is too great or too small. */
7822 /* COMMONS USED : */
7823 /* ---------------- */
7824 /* MMAPGSS,MMAPGS0,MMAPGS1,MMAPGS2. */
7825 /* ***********************************************************************
7827 /* Parameter adjustments */
7828 cgauss_dim1 = *nbpnts / 2 + 1;
7831 ibb = AdvApp2Var_SysBase::mnfndeb_();
7833 AdvApp2Var_SysBase::mgenmsg_("MMAPPTT", 7L);
7837 /* ------------------- Tests on the validity of inputs ----------------
7840 infdg = (*jordre + 1) << 1;
7841 if (*nbpnts != 8 && *nbpnts != 10 && *nbpnts != 15 && *nbpnts != 20 && *
7842 nbpnts != 25 && *nbpnts != 30 && *nbpnts != 40 && *nbpnts != 50 &&
7847 if (*jordre < -1 || *jordre > 2) {
7851 if (*ndgjac >= *nbpnts || *ndgjac < infdg) {
7855 /* --------------- Calculation of the start pointer following NBPNTS -----------
7860 iptdb += (8 - infdg) << 2;
7863 iptdb += (10 - infdg) * 5;
7866 iptdb += (15 - infdg) * 7;
7869 iptdb += (20 - infdg) * 10;
7872 iptdb += (25 - infdg) * 12;
7875 iptdb += (30 - infdg) * 15;
7878 iptdb += (40 - infdg) * 20;
7881 iptdb += (50 - infdg) * 25;
7886 ipdb0 = ipdb0 + (14 - infdg) / 2 + 1;
7889 ipdb0 = ipdb0 + (24 - infdg) / 2 + 1;
7892 /* ------------------ Choice of the common depending on JORDRE -------------
7895 if (*jordre == -1) {
7908 /* ---------------- Common MMAPGSS (case without constraints) ----------------
7912 ilong = *nbpnts / 2 << 3;
7914 for (kjac = 0; kjac <= i__1; ++kjac) {
7915 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7916 AdvApp2Var_SysBase::mcrfill_(&ilong,
7917 &mmapgss_.gslxjs[iptt - 1],
7918 &cgauss[kjac * cgauss_dim1 + 1]);
7921 /* --> Case when the number of points is uneven. */
7922 if (*nbpnts % 2 == 1) {
7925 for (kjac = 0; kjac <= i__1; kjac += 2) {
7926 cgauss[kjac * cgauss_dim1] = mmapgss_.gsl0js[iptt - 1];
7931 for (kjac = 1; kjac <= i__1; kjac += 2) {
7932 cgauss[kjac * cgauss_dim1] = 0.;
7938 /* ---------------- Common MMAPGS0 (case with constraints C0) -------------
7942 mxjac = *ndgjac - infdg;
7943 ilong = *nbpnts / 2 << 3;
7945 for (kjac = 0; kjac <= i__1; ++kjac) {
7946 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7947 AdvApp2Var_SysBase::mcrfill_(&ilong,
7948 &mmapgs0_.gslxj0[iptt - 1],
7949 &cgauss[kjac * cgauss_dim1 + 1]);
7952 /* --> Case when the number of points is uneven. */
7953 if (*nbpnts % 2 == 1) {
7956 for (kjac = 0; kjac <= i__1; kjac += 2) {
7957 cgauss[kjac * cgauss_dim1] = mmapgs0_.gsl0j0[iptt - 1];
7962 for (kjac = 1; kjac <= i__1; kjac += 2) {
7963 cgauss[kjac * cgauss_dim1] = 0.;
7969 /* ---------------- Common MMAPGS1 (case with constraints C1) -------------
7973 mxjac = *ndgjac - infdg;
7974 ilong = *nbpnts / 2 << 3;
7976 for (kjac = 0; kjac <= i__1; ++kjac) {
7977 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7978 AdvApp2Var_SysBase::mcrfill_(&ilong,
7979 &mmapgs1_.gslxj1[iptt - 1],
7980 &cgauss[kjac * cgauss_dim1 + 1]);
7983 /* --> Case when the number of points is uneven. */
7984 if (*nbpnts % 2 == 1) {
7987 for (kjac = 0; kjac <= i__1; kjac += 2) {
7988 cgauss[kjac * cgauss_dim1] = mmapgs1_.gsl0j1[iptt - 1];
7993 for (kjac = 1; kjac <= i__1; kjac += 2) {
7994 cgauss[kjac * cgauss_dim1] = 0.;
8000 /* ---------------- Common MMAPGS2 (case with constraints C2) -------------
8004 mxjac = *ndgjac - infdg;
8005 ilong = *nbpnts / 2 << 3;
8007 for (kjac = 0; kjac <= i__1; ++kjac) {
8008 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
8009 AdvApp2Var_SysBase::mcrfill_(&ilong,
8010 &mmapgs2_.gslxj2[iptt - 1],
8011 &cgauss[kjac * cgauss_dim1 + 1]);
8014 /* --> Cas of uneven number of points. */
8015 if (*nbpnts % 2 == 1) {
8018 for (kjac = 0; kjac <= i__1; kjac += 2) {
8019 cgauss[kjac * cgauss_dim1] = mmapgs2_.gsl0j2[iptt - 1];
8024 for (kjac = 1; kjac <= i__1; kjac += 2) {
8025 cgauss[kjac * cgauss_dim1] = 0.;
8031 /* ------------------------- Return the error code --------------
8033 /* --> NBPNTS is not OK */
8037 /* --> JORDRE is not OK */
8041 /* --> NDGJAC is not OK */
8046 /* -------------------------------- The end -----------------------------
8051 AdvApp2Var_SysBase::maermsg_("MMAPPTT", iercod, 7L);
8054 AdvApp2Var_SysBase::mgsomsg_("MMAPPTT", 7L);
8060 //=======================================================================
8061 //function : mmjacpt_
8063 //=======================================================================
8064 int mmjacpt_(const integer *ndimen,
8065 const integer *ncoefu,
8066 const integer *ncoefv,
8067 const integer *iordru,
8068 const integer *iordrv,
8069 const doublereal *ptclgd,
8073 /* System generated locals */
8074 integer ptccan_dim1, ptccan_dim2, ptccan_offset, ptclgd_dim1, ptclgd_dim2,
8075 ptclgd_offset, ptcaux_dim1, ptcaux_dim2, ptcaux_dim3,
8076 ptcaux_offset, i__1, i__2, i__3;
8078 /* Local variables */
8079 integer kdim, nd, ii, jj, ibb;
8081 /* ***********************************************************************
8086 /* Passage from canonical to Jacobi base for a */
8087 /* "square" in a space of arbitrary dimension. */
8091 /* SMOOTHING,BASE,LEGENDRE */
8094 /* INPUT ARGUMENTS : */
8095 /* ------------------ */
8096 /* NDIMEN : Dimension of the space. */
8097 /* NCOEFU : Degree+1 by U. */
8098 /* NCOEFV : Degree+1 by V. */
8099 /* IORDRU : Order of Jacobi polynoms by U. */
8100 /* IORDRV : Order of Jacobi polynoms by V. */
8101 /* PTCLGD : The square in the Jacobi base. */
8103 /* OUTPUT ARGUMENTS : */
8104 /* ------------------- */
8105 /* PTCAUX : Auxilliary space. */
8106 /* PTCCAN : The square in the canonic base (-1,1) */
8108 /* COMMONS USED : */
8109 /* ---------------- */
8111 /* APPLIED REFERENCES : */
8112 /* ----------------------- */
8114 /* DESCRIPTION/NOTES/LIMITATIONS : */
8115 /* ----------------------------------- */
8116 /* Cancels and replaces MJACPC */
8118 /* *********************************************************************
8120 /* Name of the routine */
8123 /* Parameter adjustments */
8124 ptccan_dim1 = *ncoefu;
8125 ptccan_dim2 = *ncoefv;
8126 ptccan_offset = ptccan_dim1 * (ptccan_dim2 + 1) + 1;
8127 ptccan -= ptccan_offset;
8128 ptcaux_dim1 = *ncoefv;
8129 ptcaux_dim2 = *ncoefu;
8130 ptcaux_dim3 = *ndimen;
8131 ptcaux_offset = ptcaux_dim1 * (ptcaux_dim2 * (ptcaux_dim3 + 1) + 1) + 1;
8132 ptcaux -= ptcaux_offset;
8133 ptclgd_dim1 = *ncoefu;
8134 ptclgd_dim2 = *ncoefv;
8135 ptclgd_offset = ptclgd_dim1 * (ptclgd_dim2 + 1) + 1;
8136 ptclgd -= ptclgd_offset;
8139 ibb = AdvApp2Var_SysBase::mnfndeb_();
8141 AdvApp2Var_SysBase::mgenmsg_("MMJACPT", 7L);
8144 /* Passage into canonical by u. */
8146 kdim = *ndimen * *ncoefv;
8147 AdvApp2Var_MathBase::mmjaccv_(ncoefu,
8150 &ptclgd[ptclgd_offset],
8151 &ptcaux[ptcaux_offset],
8152 &ptccan[ptccan_offset]);
8154 /* Swapping of u and v. */
8157 for (nd = 1; nd <= i__1; ++nd) {
8159 for (jj = 1; jj <= i__2; ++jj) {
8161 for (ii = 1; ii <= i__3; ++ii) {
8162 ptcaux[jj + (ii + (nd + ptcaux_dim3) * ptcaux_dim2) *
8163 ptcaux_dim1] = ptccan[ii + (jj + nd * ptccan_dim2) *
8172 /* Passage into canonical by v. */
8174 kdim = *ndimen * *ncoefu;
8175 AdvApp2Var_MathBase::mmjaccv_(ncoefv,
8178 &ptcaux[((ptcaux_dim3 + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1],
8179 &ptccan[ptccan_offset],
8180 &ptcaux[(((ptcaux_dim3 << 1) + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1]);
8182 /* Swapping of u and v. */
8185 for (nd = 1; nd <= i__1; ++nd) {
8187 for (jj = 1; jj <= i__2; ++jj) {
8189 for (ii = 1; ii <= i__3; ++ii) {
8190 ptccan[ii + (jj + nd * ptccan_dim2) * ptccan_dim1] = ptcaux[
8191 jj + (ii + (nd + (ptcaux_dim3 << 1)) * ptcaux_dim2) *
8200 /* ---------------------------- THAT'S ALL FOLKS ------------------------
8204 AdvApp2Var_SysBase::mgsomsg_("MMJACPT", 7L);