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 under
6 // the terms of the GNU Lesser General Public License 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;
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;
1394 ibb = AdvApp2Var_SysBase::mnfndeb_();
1396 AdvApp2Var_SysBase::mgenmsg_("MMA1JAK", 7L);
1400 /* ----------------- Recover coeffs of integration by Gauss -----------
1403 AdvApp2Var_ApproxF2var::mmapptt_(ndgjac, nbroot, iordre, cgauss, iercod);
1409 /* --------------- Calculate the curve in the base of Jacobi -----------
1412 mmmapcoe_(ndimen, ndgjac, iordre, nbroot, &somtab[somtab_offset], &diftab[
1413 diftab_offset], cgauss, &crvjac[crvjac_offset]);
1415 /* ------------------------------ The End -------------------------------
1420 AdvApp2Var_SysBase::maermsg_("MMA1JAK", iercod, 7L);
1423 AdvApp2Var_SysBase::mgsomsg_("MMA1JAK", 7L);
1428 //=======================================================================
1429 //function : mma1noc_
1431 //=======================================================================
1432 int mma1noc_(doublereal *dfuvin,
1441 /* System generated locals */
1445 /* Local variables */
1446 doublereal rider, riord;
1449 /* **********************************************************************
1454 /* Normalization of constraints of derivatives, defined on DFUVIN */
1455 /* on block DUVOUT. */
1459 /* ALL, AB_SPECIFI::VECTEUR&,DERIVEE&,NORMALISATION,&VECTEUR */
1461 /* INPUT ARGUMENTS : */
1462 /* ------------------ */
1463 /* DFUVIN: Limits of the block of definition by U and by V where
1465 /* constraints CNTRIN are defined. */
1466 /* NDIMEN: Dimension of the space. */
1467 /* IORDRE: Order of constraint imposed at the extremities of the iso. */
1468 /* (if Iso-U, it is necessary to calculate derivatives by V and vice */
1470 /* = 0, the extremities of the iso are calculated */
1471 /* = 1, additionally the 1st derivative in the direction */
1472 /* of the iso is calculated */
1473 /* = 2, additionally the 2nd derivative in the direction */
1474 /* of the iso is calculated */
1475 /* CNTRIN: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1476 /* of order IORDRE of F(Uc,v) or of F(u,Vc), following the */
1477 /* value of ISOFAV, RENORMALIZED by u and v in (-1,1). */
1478 /* DUVOUT: Limits of the block of definition by U and by V where the */
1479 /* constraints CNTOUT will be defined. */
1480 /* ISOFAV: Isoparameter fixed for the discretization; */
1481 /* = 1, discretization with fixed U=Uc and variable V. */
1482 /* = 2, discretization with fixed V=Vc and variable U. */
1483 /* IDERIV: Ordre de derivee transverse a l'iso fixee (Si Iso-U=Uc */
1484 /* is fixed, the derivative of order IDERIV is discretized by U */
1485 /* of F(Uc,v). The same if iso-V is fixed). */
1486 /* Varies from (positioning) to 2 (2nd derivative). */
1488 /* OUTPUT ARGUMENTS : */
1489 /* ------------------- */
1490 /* CNTOUT: Contains, if IORDRE>=0, IORDRE+1 derivatives */
1491 /* of order IORDRE of F(Uc,v) or of F(u,Vc), depending on the */
1492 /* value of ISOFAV, RENORMALIZED for u and v in DUVOUT. */
1494 /* COMMONS USED : */
1495 /* ---------------- */
1497 /* REFERENCES CALLED : */
1498 /* --------------------- */
1500 /* DESCRIPTION/NOTES/LIMITATIONS : */
1501 /* ------------------------------- */
1502 /* CNTRIN can be an output/input argument, */
1505 /* CALL MMA1NOC(DFUVIN,NDIMEN,IORDRE,CNTRIN,DUVOUT */
1506 /* 1 ,ISOFAV,IDERIV,CNTRIN) */
1510 /* **********************************************************************
1513 /* Name of the routine */
1516 /* Parameter adjustments */
1523 ibb = AdvApp2Var_SysBase::mnfndeb_();
1525 AdvApp2Var_SysBase::mgenmsg_("MMA1NOC", 7L);
1528 /* --------------- Determination of coefficients of normalization -------
1532 d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1533 rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1534 d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1535 riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1538 d__1 = (dfuvin[6] - dfuvin[5]) / (duvout[6] - duvout[5]);
1539 rider = AdvApp2Var_MathBase::pow__di(&d__1, ideriv);
1540 d__1 = (dfuvin[4] - dfuvin[3]) / (duvout[4] - duvout[3]);
1541 riord = AdvApp2Var_MathBase::pow__di(&d__1, iordre);
1544 /* ------------- Renormalization of the vector of constraint ---------------
1547 bid = rider * riord;
1549 for (nd = 1; nd <= i__1; ++nd) {
1550 cntout[nd] = bid * cntrin[nd];
1554 /* ------------------------------ The end -------------------------------
1558 AdvApp2Var_SysBase::mgsomsg_("MMA1NOC", 7L);
1563 //=======================================================================
1564 //function : mma1nop_
1566 //=======================================================================
1567 int mma1nop_(integer *nbroot,
1575 /* System generated locals */
1578 /* Local variables */
1579 doublereal alinu, blinu, alinv, blinv;
1582 /* ***********************************************************************
1587 /* Normalization of parameters of an iso, starting from */
1588 /* parametric block and parameters on (-1,1). */
1592 /* TOUS,AB_SPECIFI :: ISO&,POINT&,NORMALISATION,&POINT,&ISO */
1594 /* INPUT ARGUMENTS : */
1595 /* -------------------- */
1596 /* NBROOT: Nb of points of discretisation INSIDE the iso */
1597 /* defined on (-1,1). */
1598 /* ROOTLG: Table of discretization parameters on )-1,1( */
1600 /* UVFONC: Block of definition of the iso */
1601 /* ISOFAV: = 1, this is iso-u; =2, this is iso-v. */
1603 /* OUTPUT ARGUMENTS : */
1604 /* --------------------- */
1605 /* TTABLE: Table of parameters renormalized on UVFONC of the iso.
1607 /* IERCOD: = 0, OK */
1608 /* = 1, ISOFAV is out of allowed values. */
1611 /* **********************************************************************
1613 /* Name of the routine */
1616 /* Parameter adjustments */
1621 ibb = AdvApp2Var_SysBase::mnfndeb_();
1623 AdvApp2Var_SysBase::mgenmsg_("MMA1NOP", 7L);
1626 alinu = (uvfonc[4] - uvfonc[3]) / 2.;
1627 blinu = (uvfonc[4] + uvfonc[3]) / 2.;
1628 alinv = (uvfonc[6] - uvfonc[5]) / 2.;
1629 blinv = (uvfonc[6] + uvfonc[5]) / 2.;
1632 ttable[0] = uvfonc[5];
1634 for (ii = 1; ii <= i__1; ++ii) {
1635 ttable[ii] = alinv * rootlg[ii] + blinv;
1638 ttable[*nbroot + 1] = uvfonc[6];
1639 } else if (*isofav == 2) {
1640 ttable[0] = uvfonc[3];
1642 for (ii = 1; ii <= i__1; ++ii) {
1643 ttable[ii] = alinu * rootlg[ii] + blinu;
1646 ttable[*nbroot + 1] = uvfonc[4];
1653 /* ------------------------------ THE END -------------------------------
1662 AdvApp2Var_SysBase::maermsg_("MMA1NOP", iercod, 7L);
1665 AdvApp2Var_SysBase::mgsomsg_("MMA1NOP", 7L);
1672 //=======================================================================
1673 //function : mma2ac1_
1675 //=======================================================================
1676 int AdvApp2Var_ApproxF2var::mma2ac1_(integer const *ndimen,
1677 integer const *mxujac,
1678 integer const *mxvjac,
1679 integer const *iordru,
1680 integer const *iordrv,
1681 doublereal const *contr1,
1682 doublereal const * contr2,
1683 doublereal const *contr3,
1684 doublereal const *contr4,
1685 doublereal const *uhermt,
1686 doublereal const *vhermt,
1690 /* System generated locals */
1691 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
1692 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
1693 contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1,
1694 uhermt_offset, vhermt_dim1, vhermt_offset, patjac_dim1,
1695 patjac_dim2, patjac_offset, i__1, i__2, i__3, i__4, i__5;
1697 /* Local variables */
1700 doublereal bidu1, bidu2, bidv1, bidv2;
1701 integer ioru1, iorv1, ii, nd, jj, ku, kv;
1702 doublereal cnt1, cnt2, cnt3, cnt4;
1704 /* **********************************************************************
1709 /* Add polynoms of edge constraints. */
1713 /* TOUS,AB_SPECIFI::POINT&,CONTRAINTE&,ADDITION,&POLYNOME */
1715 /* INPUT ARGUMENTS : */
1716 /* ------------------ */
1717 /* NDIMEN: Dimension of the space. */
1718 /* MXUJAC: Max degree of the polynom of approximation by U. The */
1719 /* representation in the orthogonal base starts from degree */
1720 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1721 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1722 /* MXVJAC: Max degree of the polynom of approximation by V. The */
1723 /* representation in the orthogonal base starts from degree */
1724 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1725 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1726 /* IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
1727 /* to the step of constraints: C0, C1 or C2. */
1728 /* IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1729 /* to the step of constraints: C0, C1 or C2. */
1730 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
1731 /* extremities of F(U0,V0) and its derivatives. */
1732 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
1733 /* extremities of F(U1,V0) and its derivatives. */
1734 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
1735 /* extremities of F(U0,V1) and its derivatives. */
1736 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
1737 /* extremities of F(U1,V1) and its derivatives. */
1738 /* UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
1739 /* VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1740 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1741 /* of F(u,v) WITHOUT taking into account the constraints. */
1743 /* OUTPUT ARGUMENTS : */
1744 /* ------------------- */
1745 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1746 /* of F(u,v) WITH taking into account of constraints. */
1748 /* **********************************************************************
1750 /* Name of the routine */
1752 /* --------------------------- Initialization --------------------------
1755 /* Parameter adjustments */
1756 patjac_dim1 = *mxujac + 1;
1757 patjac_dim2 = *mxvjac + 1;
1758 patjac_offset = patjac_dim1 * patjac_dim2;
1759 patjac -= patjac_offset;
1760 uhermt_dim1 = (*iordru << 1) + 2;
1761 uhermt_offset = uhermt_dim1;
1762 uhermt -= uhermt_offset;
1763 vhermt_dim1 = (*iordrv << 1) + 2;
1764 vhermt_offset = vhermt_dim1;
1765 vhermt -= vhermt_offset;
1766 contr4_dim1 = *ndimen;
1767 contr4_dim2 = *iordru + 2;
1768 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
1769 contr4 -= contr4_offset;
1770 contr3_dim1 = *ndimen;
1771 contr3_dim2 = *iordru + 2;
1772 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
1773 contr3 -= contr3_offset;
1774 contr2_dim1 = *ndimen;
1775 contr2_dim2 = *iordru + 2;
1776 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
1777 contr2 -= contr2_offset;
1778 contr1_dim1 = *ndimen;
1779 contr1_dim2 = *iordru + 2;
1780 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
1781 contr1 -= contr1_offset;
1784 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1786 AdvApp2Var_SysBase::mgenmsg_("MMA2AC1", 7L);
1789 /* ------------ SUBTRACTION OF ANGULAR CONSTRAINTS -------------------
1792 ioru1 = *iordru + 1;
1793 iorv1 = *iordrv + 1;
1794 ndgu = (*iordru << 1) + 1;
1795 ndgv = (*iordrv << 1) + 1;
1798 for (jj = 1; jj <= i__1; ++jj) {
1800 for (ii = 1; ii <= i__2; ++ii) {
1802 for (nd = 1; nd <= i__3; ++nd) {
1803 cnt1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
1804 cnt2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
1805 cnt3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
1806 cnt4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
1808 for (kv = 0; kv <= i__4; ++kv) {
1809 bidv1 = vhermt[kv + ((jj << 1) - 1) * vhermt_dim1];
1810 bidv2 = vhermt[kv + (jj << 1) * vhermt_dim1];
1812 for (ku = 0; ku <= i__5; ++ku) {
1813 bidu1 = uhermt[ku + ((ii << 1) - 1) * uhermt_dim1];
1814 bidu2 = uhermt[ku + (ii << 1) * uhermt_dim1];
1815 patjac[ku + (kv + nd * patjac_dim2) * patjac_dim1] =
1816 patjac[ku + (kv + nd * patjac_dim2) *
1817 patjac_dim1] - bidu1 * bidv1 * cnt1 - bidu2 *
1818 bidv1 * cnt2 - bidu1 * bidv2 * cnt3 - bidu2 *
1831 /* ------------------------------ The end -------------------------------
1835 AdvApp2Var_SysBase::mgsomsg_("MMA2AC1", 7L);
1840 //=======================================================================
1841 //function : mma2ac2_
1843 //=======================================================================
1844 int AdvApp2Var_ApproxF2var::mma2ac2_(const integer *ndimen,
1845 const integer *mxujac,
1846 const integer *mxvjac,
1847 const integer *iordrv,
1848 const integer *nclimu,
1849 const integer *ncfiv1,
1850 const doublereal *crbiv1,
1851 const integer *ncfiv2,
1852 const doublereal *crbiv2,
1853 const doublereal *vhermt,
1857 /* System generated locals */
1858 integer crbiv1_dim1, crbiv1_dim2, crbiv1_offset, crbiv2_dim1, crbiv2_dim2,
1859 crbiv2_offset, patjac_dim1, patjac_dim2, patjac_offset,
1860 vhermt_dim1, vhermt_offset, i__1, i__2, i__3, i__4;
1862 /* Local variables */
1864 integer ndgv1, ndgv2, ii, jj, nd, kk;
1865 doublereal bid1, bid2;
1867 /* **********************************************************************
1872 /* Add polynoms of constraints */
1876 /* FUNCTION,APPROXIMATION,COEFFICIENT,POLYNOM */
1878 /* INPUT ARGUMENTS : */
1879 /* ------------------ */
1880 /* NDIMEN: Dimension of the space. */
1881 /* MXUJAC: Max degree of the polynom of approximation by U. The */
1882 /* representation in the orthogonal base starts from degree */
1883 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1884 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1885 /* MXVJAC: Max degree of the polynom of approximation by V. The */
1886 /* representation in the orthogonal base starts from degree */
1887 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
1888 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
1889 /* IORDRV: Order of the base of Jacobi (-1,0,1 or 2) by V. Corresponds */
1890 /* to the step of constraints: C0, C1 or C2. */
1891 /* NCLIMU LIMIT nb of coeff by u of the solution P(u,v)
1892 * NCFIV1: Nb of Coeff. of curves stored in CRBIV1. */
1893 /* CRBIV1: Table of coeffs of the approximation of iso-V0 and its */
1894 /* derivatives till order IORDRV. */
1895 /* NCFIV2: Nb of Coeff. of curves stored in CRBIV2. */
1896 /* CRBIV2: Table of coeffs of approximation of iso-V1 and its */
1897 /* derivatives till order IORDRV. */
1898 /* VHERMT: Coeff. of Hermit polynoms of order IORDRV. */
1899 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
1900 /* of F(u,v) WITHOUT taking into account the constraints. */
1902 /* OUTPUT ARGUMENTS : */
1903 /* ------------------- */
1904 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
1905 /* of F(u,v) WITH taking into account of constraints. */
1910 /* **********************************************************************
1912 /* Name of the routine */
1914 /* --------------------------- Initialisations --------------------------
1917 /* Parameter adjustments */
1918 patjac_dim1 = *mxujac + 1;
1919 patjac_dim2 = *mxvjac + 1;
1920 patjac_offset = patjac_dim1 * patjac_dim2;
1921 patjac -= patjac_offset;
1922 vhermt_dim1 = (*iordrv << 1) + 2;
1923 vhermt_offset = vhermt_dim1;
1924 vhermt -= vhermt_offset;
1927 crbiv2_dim1 = *nclimu;
1928 crbiv2_dim2 = *ndimen;
1929 crbiv2_offset = crbiv2_dim1 * (crbiv2_dim2 + 1);
1930 crbiv2 -= crbiv2_offset;
1931 crbiv1_dim1 = *nclimu;
1932 crbiv1_dim2 = *ndimen;
1933 crbiv1_offset = crbiv1_dim1 * (crbiv1_dim2 + 1);
1934 crbiv1 -= crbiv1_offset;
1937 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
1939 AdvApp2Var_SysBase::mgenmsg_("MMA2AC2", 7L);
1942 /* ------------ ADDING of coeff by u of curves, by v of Hermit --------
1946 for (ii = 1; ii <= i__1; ++ii) {
1947 ndgv1 = ncfiv1[ii] - 1;
1948 ndgv2 = ncfiv2[ii] - 1;
1950 for (nd = 1; nd <= i__2; ++nd) {
1951 i__3 = (*iordrv << 1) + 1;
1952 for (jj = 0; jj <= i__3; ++jj) {
1953 bid1 = vhermt[jj + ((ii << 1) - 1) * vhermt_dim1];
1955 for (kk = 0; kk <= i__4; ++kk) {
1956 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
1957 bid1 * crbiv1[kk + (nd + ii * crbiv1_dim2) *
1961 bid2 = vhermt[jj + (ii << 1) * vhermt_dim1];
1963 for (kk = 0; kk <= i__4; ++kk) {
1964 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
1965 bid2 * crbiv2[kk + (nd + ii * crbiv2_dim2) *
1976 /* ------------------------------ The end -------------------------------
1980 AdvApp2Var_SysBase::mgsomsg_("MMA2AC2", 7L);
1986 //=======================================================================
1987 //function : mma2ac3_
1989 //=======================================================================
1990 int AdvApp2Var_ApproxF2var::mma2ac3_(const integer *ndimen,
1991 const integer *mxujac,
1992 const integer *mxvjac,
1993 const integer *iordru,
1994 const integer *nclimv,
1995 const integer *ncfiu1,
1996 const doublereal * crbiu1,
1997 const integer *ncfiu2,
1998 const doublereal *crbiu2,
1999 const doublereal *uhermt,
2003 /* System generated locals */
2004 integer crbiu1_dim1, crbiu1_dim2, crbiu1_offset, crbiu2_dim1, crbiu2_dim2,
2005 crbiu2_offset, patjac_dim1, patjac_dim2, patjac_offset,
2006 uhermt_dim1, uhermt_offset, i__1, i__2, i__3, i__4;
2008 /* Local variables */
2010 integer ndgu1, ndgu2, ii, jj, nd, kk;
2011 doublereal bid1, bid2;
2013 /* **********************************************************************
2018 /* Ajout des polynomes de contraintes */
2022 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
2024 /* INPUT ARGUMENTS : */
2025 /* ------------------ */
2026 /* NDIMEN: Dimension of the space. */
2027 /* MXUJAC: Max degree of the polynom of approximation by U. The */
2028 /* representation in the orthogonal base starts from degree */
2029 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2030 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2031 /* MXVJAC: Max degree of the polynom of approximation by V. The */
2032 /* representation in the orthogonal base starts from degree */
2033 /* 0 to degree MXUJAC-2*(IORDRU+1). The polynomial base is the */
2034 /* base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
2035 /* IORDRU: Order of the base of Jacobi (-1,0,1 or 2) by U. Corresponds */
2036 /* to the step of constraints: C0, C1 or C2. */
2037 /* NCLIMV LIMIT nb of coeff by v of the solution P(u,v)
2038 * NCFIU1: Nb of Coeff. of curves stored in CRBIU1. */
2039 /* CRBIU1: Table of coeffs of the approximation of iso-U0 and its */
2040 /* derivatives till order IORDRU. */
2041 /* NCFIU2: Nb of Coeff. of curves stored in CRBIU2. */
2042 /* CRBIU2: Table of coeffs of approximation of iso-U1 and its */
2043 /* derivatives till order IORDRU */
2044 /* UHERMT: Coeff. of Hermit polynoms of order IORDRU. */
2045 /* PATJAC: Table of coefficients of the polynom P(u,v) of approximation */
2046 /* of F(u,v) WITHOUT taking into account the constraints. */
2048 /* OUTPUT ARGUMENTS : */
2049 /* ------------------- */
2050 /* PATJAC: Table of coefficients of the polynom P(u,v) by approximation */
2051 /* of F(u,v) WITH taking into account of constraints. */
2055 /* **********************************************************************
2057 /* The name of the routine */
2059 /* --------------------------- Initializations --------------------------
2062 /* Parameter adjustments */
2063 patjac_dim1 = *mxujac + 1;
2064 patjac_dim2 = *mxvjac + 1;
2065 patjac_offset = patjac_dim1 * patjac_dim2;
2066 patjac -= patjac_offset;
2067 uhermt_dim1 = (*iordru << 1) + 2;
2068 uhermt_offset = uhermt_dim1;
2069 uhermt -= uhermt_offset;
2072 crbiu2_dim1 = *nclimv;
2073 crbiu2_dim2 = *ndimen;
2074 crbiu2_offset = crbiu2_dim1 * (crbiu2_dim2 + 1);
2075 crbiu2 -= crbiu2_offset;
2076 crbiu1_dim1 = *nclimv;
2077 crbiu1_dim2 = *ndimen;
2078 crbiu1_offset = crbiu1_dim1 * (crbiu1_dim2 + 1);
2079 crbiu1 -= crbiu1_offset;
2082 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
2084 AdvApp2Var_SysBase::mgenmsg_("MMA2AC3", 7L);
2087 /* ------------ ADDING of coeff by u of curves, by v of Hermit --------
2091 for (ii = 1; ii <= i__1; ++ii) {
2092 ndgu1 = ncfiu1[ii] - 1;
2093 ndgu2 = ncfiu2[ii] - 1;
2095 for (nd = 1; nd <= i__2; ++nd) {
2097 for (jj = 0; jj <= i__3; ++jj) {
2098 bid1 = crbiu1[jj + (nd + ii * crbiu1_dim2) * crbiu1_dim1];
2099 i__4 = (*iordru << 1) + 1;
2100 for (kk = 0; kk <= i__4; ++kk) {
2101 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
2102 bid1 * uhermt[kk + ((ii << 1) - 1) * uhermt_dim1];
2108 for (jj = 0; jj <= i__3; ++jj) {
2109 bid2 = crbiu2[jj + (nd + ii * crbiu2_dim2) * crbiu2_dim1];
2110 i__4 = (*iordru << 1) + 1;
2111 for (kk = 0; kk <= i__4; ++kk) {
2112 patjac[kk + (jj + nd * patjac_dim2) * patjac_dim1] +=
2113 bid2 * uhermt[kk + (ii << 1) * uhermt_dim1];
2124 /* ------------------------------ The end -------------------------------
2128 AdvApp2Var_SysBase::mgsomsg_("MMA2AC3", 7L);
2133 //=======================================================================
2134 //function : mma2can_
2136 //=======================================================================
2137 int AdvApp2Var_ApproxF2var::mma2can_(const integer *ncfmxu,
2138 const integer *ncfmxv,
2139 const integer *ndimen,
2140 const integer *iordru,
2141 const integer *iordrv,
2142 const integer *ncoefu,
2143 const integer *ncoefv,
2144 const doublereal *patjac,
2150 /* System generated locals */
2151 integer patjac_dim1, patjac_dim2, patjac_offset, patcan_dim1, patcan_dim2,
2152 patcan_offset, i__1, i__2;
2154 /* Local variables */
2156 integer ilon1, ilon2, ii, nd;
2158 /* **********************************************************************
2163 /* Change of Jacobi base to canonical (-1,1) and writing in a greater */
2168 /* ALL,AB_SPECIFI,CARREAU&,CONVERSION,JACOBI,CANNONIQUE,&CARREAU */
2170 /* INPUT ARGUMENTS : */
2171 /* -------------------- */
2172 /* NCFMXU: Dimension by U of resulting table PATCAN */
2173 /* NCFMXV: Dimension by V of resulting table PATCAN */
2174 /* NDIMEN: Dimension of the workspace. */
2175 /* IORDRU: Order of constraint by U */
2176 /* IORDRV: Order of constraint by V. */
2177 /* NCOEFU: Nb of coeff by U of square PATJAC */
2178 /* NCOEFV: Nb of coeff by V of square PATJAC */
2179 /* PATJAC: Square in the base of Jacobi of order IORDRU by U and */
2182 /* OUTPUT ARGUMENTS : */
2183 /* --------------------- */
2184 /* PATAUX: Auxiliary Table. */
2185 /* PATCAN: Table of coefficients in the canonic base. */
2186 /* IERCOD: Error code. */
2187 /* = 0, everything goes well, and all things are equal. */
2188 /* = 1, the program refuses to process with incorrect input arguments */
2191 /* COMMONS USED : */
2192 /* ------------------ */
2194 /* REFERENCES CALLED : */
2195 /* --------------------- */
2197 /* DESCRIPTION/NOTES/LIMITATIONS : */
2198 /* ----------------------------------- */
2200 /* **********************************************************************
2204 /* Parameter adjustments */
2205 patcan_dim1 = *ncfmxu;
2206 patcan_dim2 = *ncfmxv;
2207 patcan_offset = patcan_dim1 * (patcan_dim2 + 1) + 1;
2208 patcan -= patcan_offset;
2210 patjac_dim1 = *ncoefu;
2211 patjac_dim2 = *ncoefv;
2212 patjac_offset = patjac_dim1 * (patjac_dim2 + 1) + 1;
2213 patjac -= patjac_offset;
2216 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
2218 AdvApp2Var_SysBase::mgenmsg_("MMA2CAN", 7L);
2222 if (*iordru < -1 || *iordru > 2) {
2225 if (*iordrv < -1 || *iordrv > 2) {
2228 if (*ncoefu > *ncfmxu || *ncoefv > *ncfmxv) {
2232 /* --> Pass to canonic base (-1,1) */
2233 mmjacpt_(ndimen, ncoefu, ncoefv, iordru, iordrv, &patjac[patjac_offset], &
2234 pataux[1], &patcan[patcan_offset]);
2236 /* --> Write all in a greater table */
2237 AdvApp2Var_MathBase::mmfmca8_(ncoefu,
2243 &patcan[patcan_offset],
2244 &patcan[patcan_offset]);
2246 /* --> Complete with zeros the resulting table. */
2247 ilon1 = *ncfmxu - *ncoefu;
2248 ilon2 = *ncfmxu * (*ncfmxv - *ncoefv);
2250 for (nd = 1; nd <= i__1; ++nd) {
2253 for (ii = 1; ii <= i__2; ++ii) {
2254 AdvApp2Var_SysBase::mvriraz_(&ilon1,
2255 &patcan[*ncoefu + 1 + (ii + nd * patcan_dim2) * patcan_dim1]);
2260 AdvApp2Var_SysBase::mvriraz_(&ilon2,
2261 &patcan[(*ncoefv + 1 + nd * patcan_dim2) * patcan_dim1 + 1]);
2268 /* ----------------------
2276 AdvApp2Var_SysBase::maermsg_("MMA2CAN", iercod, 7L);
2278 AdvApp2Var_SysBase::mgsomsg_("MMA2CAN", 7L);
2283 //=======================================================================
2284 //function : mma2cd1_
2286 //=======================================================================
2287 int mma2cd1_(integer *ndimen,
2310 /* System generated locals */
2311 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
2312 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
2313 contr4_dim1, contr4_dim2, contr4_offset, uhermt_dim1,
2314 uhermt_offset, vhermt_dim1, vhermt_offset, fpntbu_dim1,
2315 fpntbu_offset, fpntbv_dim1, fpntbv_offset, sosotb_dim1,
2316 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2317 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2318 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4,
2321 /* Local variables */
2322 integer ncfhu, ncfhv, nuroo, nvroo, nd, ii, jj, kk, ll, ibb, kkm,
2324 doublereal bid1, bid2, bid3, bid4;
2325 doublereal diu1, diu2, div1, div2, sou1, sou2, sov1, sov2;
2327 /* **********************************************************************
2332 /* Discretisation on the parameters of polynoms of interpolation */
2333 /* of constraints at the corners of order IORDRE. */
2337 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2339 /* INPUT ARGUMENTS : */
2340 /* ------------------ */
2341 /* NDIMEN: Dimension of the space. */
2342 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2343 /* This is also the nb of root of Legendre polynom where discretization is done. */
2344 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
2346 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
2347 /* This is also the nb of root of Legendre polynom where discretization is done. */
2348 /* VROOTL: Table of discretization parameters on (-1,1) by V. */
2349 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
2350 /* = 0, calculate the extremities of iso-V */
2351 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2352 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2353 /* IORDRV: Order of constraint imposed at the extremities of iso-U */
2354 /* = 0, calculate the extremities of iso-U */
2355 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
2356 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
2357 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
2358 /* extremities of F(U0,V0) and its derivatives. */
2359 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
2360 /* extremities of F(U1,V0) and its derivatives. */
2361 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
2362 /* extremities of F(U0,V1) and its derivatives. */
2363 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
2364 /* extremities of F(U1,V1) and its derivatives. */
2365 /* SOSOTB: Preinitialized table (input/output argument). */
2366 /* DISOTB: Preinitialized table (input/output argument). */
2367 /* SODITB: Preinitialized table (input/output argument). */
2368 /* DIDITB: Preinitialized table (input/output argument) */
2370 /* OUTPUT ARGUMENTS : */
2371 /* ------------------- */
2372 /* FPNTBU: Auxiliary table. */
2373 /* FPNTBV: Auxiliary table. */
2374 /* UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
2375 /* VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2376 /* SOSOTB: Table where the terms of constraints are added */
2377 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2378 /* with ui and vj positive roots of the Legendre polynom */
2379 /* of degree NBPNTU and NBPNTV respectively. */
2380 /* DISOTB: Table where the terms of constraints are added */
2381 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2382 /* with ui and vj positive roots of the polynom of Legendre */
2383 /* of degree NBPNTU and NBPNTV respectively. */
2384 /* SODITB: Table where the terms of constraints are added */
2385 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2386 /* with ui and vj positive roots of the polynom of Legendre */
2387 /* of degree NBPNTU and NBPNTV respectively. */
2388 /* DIDITB: Table where the terms of constraints are added */
2389 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2390 /* with ui and vj positive roots of the polynom of Legendre */
2391 /* of degree NBPNTU and NBPNTV respectively. */
2393 /* COMMONS USED : */
2394 /* ---------------- */
2396 /* REFERENCES CALLED : */
2397 /* ----------------------- */
2399 /* DESCRIPTION/NOTES/LIMITATIONS : */
2400 /* ----------------------------------- */
2403 /* **********************************************************************
2406 /* Name of the routine */
2409 /* Parameter adjustments */
2411 diditb_dim1 = *nbpntu / 2 + 1;
2412 diditb_dim2 = *nbpntv / 2 + 1;
2413 diditb_offset = diditb_dim1 * diditb_dim2;
2414 diditb -= diditb_offset;
2415 disotb_dim1 = *nbpntu / 2;
2416 disotb_dim2 = *nbpntv / 2;
2417 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2418 disotb -= disotb_offset;
2419 soditb_dim1 = *nbpntu / 2;
2420 soditb_dim2 = *nbpntv / 2;
2421 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2422 soditb -= soditb_offset;
2423 sosotb_dim1 = *nbpntu / 2 + 1;
2424 sosotb_dim2 = *nbpntv / 2 + 1;
2425 sosotb_offset = sosotb_dim1 * sosotb_dim2;
2426 sosotb -= sosotb_offset;
2428 uhermt_dim1 = (*iordru << 1) + 2;
2429 uhermt_offset = uhermt_dim1;
2430 uhermt -= uhermt_offset;
2431 fpntbu_dim1 = *nbpntu;
2432 fpntbu_offset = fpntbu_dim1 + 1;
2433 fpntbu -= fpntbu_offset;
2434 vhermt_dim1 = (*iordrv << 1) + 2;
2435 vhermt_offset = vhermt_dim1;
2436 vhermt -= vhermt_offset;
2437 fpntbv_dim1 = *nbpntv;
2438 fpntbv_offset = fpntbv_dim1 + 1;
2439 fpntbv -= fpntbv_offset;
2440 contr4_dim1 = *ndimen;
2441 contr4_dim2 = *iordru + 2;
2442 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
2443 contr4 -= contr4_offset;
2444 contr3_dim1 = *ndimen;
2445 contr3_dim2 = *iordru + 2;
2446 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
2447 contr3 -= contr3_offset;
2448 contr2_dim1 = *ndimen;
2449 contr2_dim2 = *iordru + 2;
2450 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
2451 contr2 -= contr2_offset;
2452 contr1_dim1 = *ndimen;
2453 contr1_dim2 = *iordru + 2;
2454 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
2455 contr1 -= contr1_offset;
2458 ibb = AdvApp2Var_SysBase::mnfndeb_();
2460 AdvApp2Var_SysBase::mgenmsg_("MMA2CD1", 7L);
2463 /* ------------------- Discretisation of Hermite polynoms -----------
2466 ncfhu = (*iordru + 1) << 1;
2468 for (ii = 1; ii <= i__1; ++ii) {
2470 for (ll = 1; ll <= i__2; ++ll) {
2471 AdvApp2Var_MathBase::mmmpocur_(&ncfhu, &c__1, &ncfhu, &uhermt[ii * uhermt_dim1], &
2472 urootl[ll], &fpntbu[ll + ii * fpntbu_dim1]);
2477 ncfhv = (*iordrv + 1) << 1;
2479 for (jj = 1; jj <= i__1; ++jj) {
2481 for (kk = 1; kk <= i__2; ++kk) {
2482 AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[jj * vhermt_dim1], &
2483 vrootl[kk], &fpntbv[kk + jj * fpntbv_dim1]);
2489 /* ---- The discretizations of polynoms of constraints are subtracted ----
2492 nuroo = *nbpntu / 2;
2493 nvroo = *nbpntv / 2;
2495 for (nd = 1; nd <= i__1; ++nd) {
2498 for (jj = 1; jj <= i__2; ++jj) {
2500 for (ii = 1; ii <= i__3; ++ii) {
2501 bid1 = contr1[nd + (ii + jj * contr1_dim2) * contr1_dim1];
2502 bid2 = contr2[nd + (ii + jj * contr2_dim2) * contr2_dim1];
2503 bid3 = contr3[nd + (ii + jj * contr3_dim2) * contr3_dim1];
2504 bid4 = contr4[nd + (ii + jj * contr4_dim2) * contr4_dim1];
2507 for (kk = 1; kk <= i__4; ++kk) {
2508 kkp = (*nbpntv + 1) / 2 + kk;
2509 kkm = nvroo - kk + 1;
2510 sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] +
2511 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2512 div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] -
2513 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2514 sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[kkm
2515 + (jj << 1) * fpntbv_dim1];
2516 div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[kkm
2517 + (jj << 1) * fpntbv_dim1];
2519 for (ll = 1; ll <= i__5; ++ll) {
2520 llp = (*nbpntu + 1) / 2 + ll;
2521 llm = nuroo - ll + 1;
2522 sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] +
2523 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2524 diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] -
2525 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2526 sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2527 llm + (ii << 1) * fpntbu_dim1];
2528 diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2529 llm + (ii << 1) * fpntbu_dim1];
2530 sosotb[ll + (kk + nd * sosotb_dim2) * sosotb_dim1] =
2531 sosotb[ll + (kk + nd * sosotb_dim2) *
2532 sosotb_dim1] - bid1 * sou1 * sov1 - bid2 *
2533 sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2535 soditb[ll + (kk + nd * soditb_dim2) * soditb_dim1] =
2536 soditb[ll + (kk + nd * soditb_dim2) *
2537 soditb_dim1] - bid1 * sou1 * div1 - bid2 *
2538 sou2 * div1 - bid3 * sou1 * div2 - bid4 *
2540 disotb[ll + (kk + nd * disotb_dim2) * disotb_dim1] =
2541 disotb[ll + (kk + nd * disotb_dim2) *
2542 disotb_dim1] - bid1 * diu1 * sov1 - bid2 *
2543 diu2 * sov1 - bid3 * diu1 * sov2 - bid4 *
2545 diditb[ll + (kk + nd * diditb_dim2) * diditb_dim1] =
2546 diditb[ll + (kk + nd * diditb_dim2) *
2547 diditb_dim1] - bid1 * diu1 * div1 - bid2 *
2548 diu2 * div1 - bid3 * diu1 * div2 - bid4 *
2555 /* ------------ Case when the discretization is done only on the roots
2557 /* ---------- of Legendre polynom of uneven degree, 0 is root
2560 if (*nbpntu % 2 == 1) {
2561 sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2562 sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2564 for (kk = 1; kk <= i__4; ++kk) {
2565 kkp = (*nbpntv + 1) / 2 + kk;
2566 kkm = nvroo - kk + 1;
2567 sov1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] +
2568 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2569 div1 = fpntbv[kkp + ((jj << 1) - 1) * fpntbv_dim1] -
2570 fpntbv[kkm + ((jj << 1) - 1) * fpntbv_dim1];
2571 sov2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] + fpntbv[
2572 kkm + (jj << 1) * fpntbv_dim1];
2573 div2 = fpntbv[kkp + (jj << 1) * fpntbv_dim1] - fpntbv[
2574 kkm + (jj << 1) * fpntbv_dim1];
2575 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1] =
2576 sosotb[(kk + nd * sosotb_dim2) * sosotb_dim1]
2577 - bid1 * sou1 * sov1 - bid2 * sou2 * sov1 -
2578 bid3 * sou1 * sov2 - bid4 * sou2 * sov2;
2579 diditb[(kk + nd * diditb_dim2) * diditb_dim1] =
2580 diditb[(kk + nd * diditb_dim2) * diditb_dim1]
2581 - bid1 * sou1 * div1 - bid2 * sou2 * div1 -
2582 bid3 * sou1 * div2 - bid4 * sou2 * div2;
2587 if (*nbpntv % 2 == 1) {
2588 sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2589 sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2591 for (ll = 1; ll <= i__4; ++ll) {
2592 llp = (*nbpntu + 1) / 2 + ll;
2593 llm = nuroo - ll + 1;
2594 sou1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] +
2595 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2596 diu1 = fpntbu[llp + ((ii << 1) - 1) * fpntbu_dim1] -
2597 fpntbu[llm + ((ii << 1) - 1) * fpntbu_dim1];
2598 sou2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] + fpntbu[
2599 llm + (ii << 1) * fpntbu_dim1];
2600 diu2 = fpntbu[llp + (ii << 1) * fpntbu_dim1] - fpntbu[
2601 llm + (ii << 1) * fpntbu_dim1];
2602 sosotb[ll + nd * sosotb_dim2 * sosotb_dim1] = sosotb[
2603 ll + nd * sosotb_dim2 * sosotb_dim1] - bid1 *
2604 sou1 * sov1 - bid2 * sou2 * sov1 - bid3 *
2605 sou1 * sov2 - bid4 * sou2 * sov2;
2606 diditb[ll + nd * diditb_dim2 * diditb_dim1] = diditb[
2607 ll + nd * diditb_dim2 * diditb_dim1] - bid1 *
2608 diu1 * sov1 - bid2 * diu2 * sov1 - bid3 *
2609 diu1 * sov2 - bid4 * diu2 * sov2;
2614 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2615 sou1 = fpntbu[nuroo + 1 + ((ii << 1) - 1) * fpntbu_dim1];
2616 sou2 = fpntbu[nuroo + 1 + (ii << 1) * fpntbu_dim1];
2617 sov1 = fpntbv[nvroo + 1 + ((jj << 1) - 1) * fpntbv_dim1];
2618 sov2 = fpntbv[nvroo + 1 + (jj << 1) * fpntbv_dim1];
2619 sosotb[nd * sosotb_dim2 * sosotb_dim1] = sosotb[nd *
2620 sosotb_dim2 * sosotb_dim1] - bid1 * sou1 * sov1 -
2621 bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2623 diditb[nd * diditb_dim2 * diditb_dim1] = diditb[nd *
2624 diditb_dim2 * diditb_dim1] - bid1 * sou1 * sov1 -
2625 bid2 * sou2 * sov1 - bid3 * sou1 * sov2 - bid4 *
2637 /* ------------------------------ The End -------------------------------
2642 AdvApp2Var_SysBase::mgsomsg_("MMA2CD1", 7L);
2647 //=======================================================================
2648 //function : mma2cd2_
2650 //=======================================================================
2651 int mma2cd2_(integer *ndimen,
2669 /* System generated locals */
2670 integer sotbv1_dim1, sotbv1_dim2, sotbv1_offset, sotbv2_dim1, sotbv2_dim2,
2671 sotbv2_offset, ditbv1_dim1, ditbv1_dim2, ditbv1_offset,
2672 ditbv2_dim1, ditbv2_dim2, ditbv2_offset, fpntab_dim1,
2673 fpntab_offset, vhermt_dim1, vhermt_offset, sosotb_dim1,
2674 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2675 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2676 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2678 /* Local variables */
2679 integer ncfhv, nuroo, nvroo, ii, nd, jj, kk, ibb, jjm, jjp;
2680 doublereal bid1, bid2, bid3, bid4;
2682 /* **********************************************************************
2686 /* Discretisation on the parameters of polynoms of interpolation */
2687 /* of constraints on 2 borders iso-V of order IORDRV. */
2692 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
2696 /* INPUT ARGUMENTS : */
2697 /* ------------------ */
2698 /* NDIMEN: Dimension of the space. */
2699 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
2700 /* This is also the nb of root of Legendre polynom where discretization is done. */
2701 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
2703 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
2704 /* This is also the nb of root of Legendre polynom where discretization is done. */
2705 /* VROOTL: Table of discretization parameters on (-1,1) by V. */
2706 /* IORDRV: Order of constraint imposed at the extremities of iso-V */
2707 /* = 0, calculate the extremities of iso-V */
2708 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
2709 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
2710 /* SOTBV1: Table of NBPNTV/2 sums of 2 index points */
2711 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2712 /* SOTBV2: Table of NBPNTV/2 sums of 2 index points */
2713 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2714 /* DITBV1: Table of NBPNTV/2 differences of 2 index points */
2715 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
2716 /* DITBV2: Table of NBPNTV/2 differences of 2 index points */
2717 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
2718 /* SOSOTB: Preinitialized table (input/output argument). */
2719 /* DISOTB: Preinitialized table (input/output argument). */
2720 /* SODITB: Preinitialized table (input/output argument). */
2721 /* DIDITB: Preinitialized table (input/output argument) */
2723 /* OUTPUT ARGUMENTS : */
2724 /* ------------------- */
2725 /* FPNTAB: Auxiliary table. */
2726 /* VHERMT: Table of 2*(IORDRV+1) coeff. of 2*(IORDRV+1) polynoms of Hermite. */
2727 /* SOSOTB: Table where the terms of constraints are added */
2728 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
2729 /* with ui and vj positive roots of the Legendre polynom */
2730 /* of degree NBPNTU and NBPNTV respectively. */
2731 /* DISOTB: Table where the terms of constraints are added */
2732 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
2733 /* with ui and vj positive roots of the polynom of Legendre */
2734 /* of degree NBPNTU and NBPNTV respectively. */
2735 /* SODITB: Table where the terms of constraints are added */
2736 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
2737 /* with ui and vj positive roots of the polynom of Legendre */
2738 /* of degree NBPNTU and NBPNTV respectively. */
2739 /* DIDITB: Table where the terms of constraints are added */
2740 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
2741 /* with ui and vj positive roots of the polynom of Legendre */
2742 /* of degree NBPNTU and NBPNTV respectively. */
2744 /* COMMONS USED : */
2745 /* ---------------- */
2747 /* REFERENCES CALLED : */
2748 /* ----------------------- */
2750 /* DESCRIPTION/NOTES/LIMITATIONS : */
2751 /* ----------------------------------- */
2755 /* **********************************************************************
2758 /* Name of the routine */
2761 /* Parameter adjustments */
2762 diditb_dim1 = *nbpntu / 2 + 1;
2763 diditb_dim2 = *nbpntv / 2 + 1;
2764 diditb_offset = diditb_dim1 * diditb_dim2;
2765 diditb -= diditb_offset;
2766 disotb_dim1 = *nbpntu / 2;
2767 disotb_dim2 = *nbpntv / 2;
2768 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
2769 disotb -= disotb_offset;
2770 soditb_dim1 = *nbpntu / 2;
2771 soditb_dim2 = *nbpntv / 2;
2772 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
2773 soditb -= soditb_offset;
2774 sosotb_dim1 = *nbpntu / 2 + 1;
2775 sosotb_dim2 = *nbpntv / 2 + 1;
2776 sosotb_offset = sosotb_dim1 * sosotb_dim2;
2777 sosotb -= sosotb_offset;
2779 vhermt_dim1 = (*iordrv << 1) + 2;
2780 vhermt_offset = vhermt_dim1;
2781 vhermt -= vhermt_offset;
2782 fpntab_dim1 = *nbpntv;
2783 fpntab_offset = fpntab_dim1 + 1;
2784 fpntab -= fpntab_offset;
2785 ditbv2_dim1 = *nbpntu / 2 + 1;
2786 ditbv2_dim2 = *ndimen;
2787 ditbv2_offset = ditbv2_dim1 * (ditbv2_dim2 + 1);
2788 ditbv2 -= ditbv2_offset;
2789 ditbv1_dim1 = *nbpntu / 2 + 1;
2790 ditbv1_dim2 = *ndimen;
2791 ditbv1_offset = ditbv1_dim1 * (ditbv1_dim2 + 1);
2792 ditbv1 -= ditbv1_offset;
2793 sotbv2_dim1 = *nbpntu / 2 + 1;
2794 sotbv2_dim2 = *ndimen;
2795 sotbv2_offset = sotbv2_dim1 * (sotbv2_dim2 + 1);
2796 sotbv2 -= sotbv2_offset;
2797 sotbv1_dim1 = *nbpntu / 2 + 1;
2798 sotbv1_dim2 = *ndimen;
2799 sotbv1_offset = sotbv1_dim1 * (sotbv1_dim2 + 1);
2800 sotbv1 -= sotbv1_offset;
2803 ibb = AdvApp2Var_SysBase::mnfndeb_();
2805 AdvApp2Var_SysBase::mgenmsg_("MMA2CD2", 7L);
2808 /* ------------------- Discretization of Hermit polynoms -----------
2811 ncfhv = (*iordrv + 1) << 1;
2813 for (ii = 1; ii <= i__1; ++ii) {
2815 for (jj = 1; jj <= i__2; ++jj) {
2816 AdvApp2Var_MathBase::mmmpocur_(&ncfhv, &c__1, &ncfhv, &vhermt[ii * vhermt_dim1], &
2817 vrootl[jj], &fpntab[jj + ii * fpntab_dim1]);
2823 /* ---- The discretizations of polynoms of constraints are subtracted ----
2826 nuroo = *nbpntu / 2;
2827 nvroo = *nbpntv / 2;
2830 for (nd = 1; nd <= i__1; ++nd) {
2832 for (ii = 1; ii <= i__2; ++ii) {
2835 for (kk = 1; kk <= i__3; ++kk) {
2836 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1];
2837 bid2 = sotbv2[kk + (nd + ii * sotbv2_dim2) * sotbv2_dim1];
2838 bid3 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1];
2839 bid4 = ditbv2[kk + (nd + ii * ditbv2_dim2) * ditbv2_dim1];
2841 for (jj = 1; jj <= i__4; ++jj) {
2842 jjp = (*nbpntv + 1) / 2 + jj;
2843 jjm = nvroo - jj + 1;
2844 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] =
2845 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
2846 - bid1 * (fpntab[jjp + ((ii << 1) - 1) *
2847 fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) *
2848 fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) *
2849 fpntab_dim1] + fpntab[jjm + (ii << 1) *
2851 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] =
2852 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
2853 - bid3 * (fpntab[jjp + ((ii << 1) - 1) *
2854 fpntab_dim1] + fpntab[jjm + ((ii << 1) - 1) *
2855 fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) *
2856 fpntab_dim1] + fpntab[jjm + (ii << 1) *
2858 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] =
2859 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
2860 - bid1 * (fpntab[jjp + ((ii << 1) - 1) *
2861 fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) *
2862 fpntab_dim1]) - bid2 * (fpntab[jjp + (ii << 1) *
2863 fpntab_dim1] - fpntab[jjm + (ii << 1) *
2865 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] =
2866 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
2867 - bid3 * (fpntab[jjp + ((ii << 1) - 1) *
2868 fpntab_dim1] - fpntab[jjm + ((ii << 1) - 1) *
2869 fpntab_dim1]) - bid4 * (fpntab[jjp + (ii << 1) *
2870 fpntab_dim1] - fpntab[jjm + (ii << 1) *
2879 /* ------------ Case when the discretization is done only on the roots */
2880 /* ---------- of Legendre polynom of uneven degree, 0 is root */
2883 if (*nbpntv % 2 == 1) {
2885 for (ii = 1; ii <= i__2; ++ii) {
2887 for (kk = 1; kk <= i__3; ++kk) {
2888 bid1 = sotbv1[kk + (nd + ii * sotbv1_dim2) * sotbv1_dim1]
2889 * fpntab[nvroo + 1 + ((ii << 1) - 1) *
2890 fpntab_dim1] + sotbv2[kk + (nd + ii * sotbv2_dim2)
2891 * sotbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) *
2893 sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2894 bid2 = ditbv1[kk + (nd + ii * ditbv1_dim2) * ditbv1_dim1]
2895 * fpntab[nvroo + 1 + ((ii << 1) - 1) *
2896 fpntab_dim1] + ditbv2[kk + (nd + ii * ditbv2_dim2)
2897 * ditbv2_dim1] * fpntab[nvroo + 1 + (ii << 1) *
2899 diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
2906 if (*nbpntu % 2 == 1) {
2908 for (ii = 1; ii <= i__2; ++ii) {
2910 for (jj = 1; jj <= i__3; ++jj) {
2911 jjp = (*nbpntv + 1) / 2 + jj;
2912 jjm = nvroo - jj + 1;
2913 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2914 fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] +
2915 fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) +
2916 sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2917 fpntab[jjp + (ii << 1) * fpntab_dim1] + fpntab[
2918 jjm + (ii << 1) * fpntab_dim1]);
2919 sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
2920 bid2 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * (
2921 fpntab[jjp + ((ii << 1) - 1) * fpntab_dim1] -
2922 fpntab[jjm + ((ii << 1) - 1) * fpntab_dim1]) +
2923 sotbv2[(nd + ii * sotbv2_dim2) * sotbv2_dim1] * (
2924 fpntab[jjp + (ii << 1) * fpntab_dim1] - fpntab[
2925 jjm + (ii << 1) * fpntab_dim1]);
2926 diditb[jj + nd * diditb_dim2 * diditb_dim1] -= bid2;
2933 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
2935 for (ii = 1; ii <= i__2; ++ii) {
2936 bid1 = sotbv1[(nd + ii * sotbv1_dim2) * sotbv1_dim1] * fpntab[
2937 nvroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbv2[(
2938 nd + ii * sotbv2_dim2) * sotbv2_dim1] * fpntab[nvroo
2939 + 1 + (ii << 1) * fpntab_dim1];
2940 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
2949 /* ------------------------------ The End -------------------------------
2954 AdvApp2Var_SysBase::mgsomsg_("MMA2CD2", 7L);
2959 //=======================================================================
2960 //function : mma2cd3_
2962 //=======================================================================
2963 int mma2cd3_(integer *ndimen,
2982 /* System generated locals */
2983 integer sotbu1_dim1, sotbu1_dim2, sotbu1_offset, sotbu2_dim1, sotbu2_dim2,
2984 sotbu2_offset, ditbu1_dim1, ditbu1_dim2, ditbu1_offset,
2985 ditbu2_dim1, ditbu2_dim2, ditbu2_offset, fpntab_dim1,
2986 fpntab_offset, uhermt_dim1, uhermt_offset, sosotb_dim1,
2987 sosotb_dim2, sosotb_offset, diditb_dim1, diditb_dim2,
2988 diditb_offset, soditb_dim1, soditb_dim2, soditb_offset,
2989 disotb_dim1, disotb_dim2, disotb_offset, i__1, i__2, i__3, i__4;
2991 /* Local variables */
2992 integer ncfhu, nuroo, nvroo, ii, nd, jj, kk, ibb, kkm, kkp;
2993 doublereal bid1, bid2, bid3, bid4;
2995 /* **********************************************************************
2999 /* Discretisation on the parameters of polynoms of interpolation */
3000 /* of constraints on 2 borders iso-U of order IORDRU. */
3005 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3007 /* INPUT ARGUMENTS : */
3008 /* ------------------ */
3009 /* NDIMEN: Dimension of the space. */
3010 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3011 /* This is also the nb of root of Legendre polynom where discretization is done. */
3012 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
3014 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
3015 /* This is also the nb of root of Legendre polynom where discretization is done. */
3016 /* IORDRV: Order of constraint imposed at the extremities of iso-V */
3017 /* = 0, calculate the extremities of iso-V */
3018 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3019 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3020 /* SOTBU1: Table of NBPNTU/2 sums of 2 index points */
3021 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3022 /* SOTBU2: Table of NBPNTV/2 sums of 2 index points */
3023 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3024 /* DITBU1: Table of NBPNTU/2 differences of 2 index points */
3025 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3026 /* DITBU2: Table of NBPNTU/2 differences of 2 index points */
3027 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3028 /* SOSOTB: Preinitialized table (input/output argument). */
3029 /* DISOTB: Preinitialized table (input/output argument). */
3030 /* SODITB: Preinitialized table (input/output argument). */
3031 /* DIDITB: Preinitialized table (input/output argument) */
3033 /* OUTPUT ARGUMENTS : */
3034 /* ------------------- */
3035 /* FPNTAB: Auxiliary table. */
3036 /* UHERMT: Table of 2*(IORDRU+1) coeff. of 2*(IORDRU+1) polynoms of Hermite. */
3037 /* SOSOTB: Table where the terms of constraints are added */
3038 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3039 /* with ui and vj positive roots of the Legendre polynom */
3040 /* of degree NBPNTU and NBPNTV respectively. */
3041 /* DISOTB: Table where the terms of constraints are added */
3042 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3043 /* with ui and vj positive roots of the polynom of Legendre */
3044 /* of degree NBPNTU and NBPNTV respectively. */
3045 /* SODITB: Table where the terms of constraints are added */
3046 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3047 /* with ui and vj positive roots of the polynom of Legendre */
3048 /* of degree NBPNTU and NBPNTV respectively. */
3049 /* DIDITB: Table where the terms of constraints are added */
3050 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3051 /* with ui and vj positive roots of the polynom of Legendre */
3052 /* of degree NBPNTU and NBPNTV respectively. */
3054 /* COMMONS USED : */
3055 /* ---------------- */
3057 /* REFERENCES CALLED : */
3058 /* ----------------------- */
3060 /* DESCRIPTION/NOTES/LIMITATIONS : */
3061 /* ----------------------------------- */
3063 /* $ HISTORIQUE DES MODIFICATIONS : */
3064 /* -------------------------------- */
3065 /* 08-08-1991: RBD; Creation. */
3067 /* **********************************************************************
3070 /* Name of the routine */
3073 /* Parameter adjustments */
3075 diditb_dim1 = *nbpntu / 2 + 1;
3076 diditb_dim2 = *nbpntv / 2 + 1;
3077 diditb_offset = diditb_dim1 * diditb_dim2;
3078 diditb -= diditb_offset;
3079 disotb_dim1 = *nbpntu / 2;
3080 disotb_dim2 = *nbpntv / 2;
3081 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3082 disotb -= disotb_offset;
3083 soditb_dim1 = *nbpntu / 2;
3084 soditb_dim2 = *nbpntv / 2;
3085 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3086 soditb -= soditb_offset;
3087 sosotb_dim1 = *nbpntu / 2 + 1;
3088 sosotb_dim2 = *nbpntv / 2 + 1;
3089 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3090 sosotb -= sosotb_offset;
3091 uhermt_dim1 = (*iordru << 1) + 2;
3092 uhermt_offset = uhermt_dim1;
3093 uhermt -= uhermt_offset;
3094 fpntab_dim1 = *nbpntu;
3095 fpntab_offset = fpntab_dim1 + 1;
3096 fpntab -= fpntab_offset;
3097 ditbu2_dim1 = *nbpntv / 2 + 1;
3098 ditbu2_dim2 = *ndimen;
3099 ditbu2_offset = ditbu2_dim1 * (ditbu2_dim2 + 1);
3100 ditbu2 -= ditbu2_offset;
3101 ditbu1_dim1 = *nbpntv / 2 + 1;
3102 ditbu1_dim2 = *ndimen;
3103 ditbu1_offset = ditbu1_dim1 * (ditbu1_dim2 + 1);
3104 ditbu1 -= ditbu1_offset;
3105 sotbu2_dim1 = *nbpntv / 2 + 1;
3106 sotbu2_dim2 = *ndimen;
3107 sotbu2_offset = sotbu2_dim1 * (sotbu2_dim2 + 1);
3108 sotbu2 -= sotbu2_offset;
3109 sotbu1_dim1 = *nbpntv / 2 + 1;
3110 sotbu1_dim2 = *ndimen;
3111 sotbu1_offset = sotbu1_dim1 * (sotbu1_dim2 + 1);
3112 sotbu1 -= sotbu1_offset;
3115 ibb = AdvApp2Var_SysBase::mnfndeb_();
3117 AdvApp2Var_SysBase::mgenmsg_("MMA2CD3", 7L);
3120 /* ------------------- Discretization of polynoms of Hermit -----------
3123 ncfhu = (*iordru + 1) << 1;
3125 for (ii = 1; ii <= i__1; ++ii) {
3127 for (kk = 1; kk <= i__2; ++kk) {
3128 AdvApp2Var_MathBase::mmmpocur_(&ncfhu,
3131 &uhermt[ii * uhermt_dim1],
3133 &fpntab[kk + ii * fpntab_dim1]);
3139 /* ---- The discretizations of polynoms of constraints are subtracted ----
3142 nvroo = *nbpntv / 2;
3143 nuroo = *nbpntu / 2;
3146 for (nd = 1; nd <= i__1; ++nd) {
3148 for (ii = 1; ii <= i__2; ++ii) {
3151 for (jj = 1; jj <= i__3; ++jj) {
3152 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1];
3153 bid2 = sotbu2[jj + (nd + ii * sotbu2_dim2) * sotbu2_dim1];
3154 bid3 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1];
3155 bid4 = ditbu2[jj + (nd + ii * ditbu2_dim2) * ditbu2_dim1];
3157 for (kk = 1; kk <= i__4; ++kk) {
3158 kkp = (*nbpntu + 1) / 2 + kk;
3159 kkm = nuroo - kk + 1;
3160 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1] =
3161 sosotb[kk + (jj + nd * sosotb_dim2) * sosotb_dim1]
3162 - bid1 * (fpntab[kkp + ((ii << 1) - 1) *
3163 fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) *
3164 fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) *
3165 fpntab_dim1] + fpntab[kkm + (ii << 1) *
3167 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1] =
3168 disotb[kk + (jj + nd * disotb_dim2) * disotb_dim1]
3169 - bid1 * (fpntab[kkp + ((ii << 1) - 1) *
3170 fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) *
3171 fpntab_dim1]) - bid2 * (fpntab[kkp + (ii << 1) *
3172 fpntab_dim1] - fpntab[kkm + (ii << 1) *
3174 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1] =
3175 soditb[kk + (jj + nd * soditb_dim2) * soditb_dim1]
3176 - bid3 * (fpntab[kkp + ((ii << 1) - 1) *
3177 fpntab_dim1] + fpntab[kkm + ((ii << 1) - 1) *
3178 fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) *
3179 fpntab_dim1] + fpntab[kkm + (ii << 1) *
3181 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1] =
3182 diditb[kk + (jj + nd * diditb_dim2) * diditb_dim1]
3183 - bid3 * (fpntab[kkp + ((ii << 1) - 1) *
3184 fpntab_dim1] - fpntab[kkm + ((ii << 1) - 1) *
3185 fpntab_dim1]) - bid4 * (fpntab[kkp + (ii << 1) *
3186 fpntab_dim1] - fpntab[kkm + (ii << 1) *
3195 /* ------------ Case when the discretization is done only on the roots */
3196 /* ---------- of Legendre polynom of uneven degree, 0 is root */
3200 if (*nbpntu % 2 == 1) {
3202 for (ii = 1; ii <= i__2; ++ii) {
3204 for (jj = 1; jj <= i__3; ++jj) {
3205 bid1 = sotbu1[jj + (nd + ii * sotbu1_dim2) * sotbu1_dim1]
3206 * fpntab[nuroo + 1 + ((ii << 1) - 1) *
3207 fpntab_dim1] + sotbu2[jj + (nd + ii * sotbu2_dim2)
3208 * sotbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) *
3210 sosotb[(jj + nd * sosotb_dim2) * sosotb_dim1] -= bid1;
3211 bid2 = ditbu1[jj + (nd + ii * ditbu1_dim2) * ditbu1_dim1]
3212 * fpntab[nuroo + 1 + ((ii << 1) - 1) *
3213 fpntab_dim1] + ditbu2[jj + (nd + ii * ditbu2_dim2)
3214 * ditbu2_dim1] * fpntab[nuroo + 1 + (ii << 1) *
3216 diditb[(jj + nd * diditb_dim2) * diditb_dim1] -= bid2;
3223 if (*nbpntv % 2 == 1) {
3225 for (ii = 1; ii <= i__2; ++ii) {
3227 for (kk = 1; kk <= i__3; ++kk) {
3228 kkp = (*nbpntu + 1) / 2 + kk;
3229 kkm = nuroo - kk + 1;
3230 bid1 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * (
3231 fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] +
3232 fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1]) +
3233 sotbu2[(nd + ii * sotbu2_dim2) * sotbu2_dim1] * (
3234 fpntab[kkp + (ii << 1) * fpntab_dim1] + fpntab[
3235 kkm + (ii << 1) * fpntab_dim1]);
3236 sosotb[kk + nd * sosotb_dim2 * sosotb_dim1] -= bid1;
3237 bid2 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * (
3238 fpntab[kkp + ((ii << 1) - 1) * fpntab_dim1] -
3239 fpntab[kkm + ((ii << 1) - 1) * fpntab_dim1]) +
3240 sotbu2[(nd + ii * sotbu2_dim2) * sotbu2_dim1] * (
3241 fpntab[kkp + (ii << 1) * fpntab_dim1] - fpntab[
3242 kkm + (ii << 1) * fpntab_dim1]);
3243 diditb[kk + nd * diditb_dim2 * diditb_dim1] -= bid2;
3250 if (*nbpntu % 2 == 1 && *nbpntv % 2 == 1) {
3252 for (ii = 1; ii <= i__2; ++ii) {
3253 bid1 = sotbu1[(nd + ii * sotbu1_dim2) * sotbu1_dim1] * fpntab[
3254 nuroo + 1 + ((ii << 1) - 1) * fpntab_dim1] + sotbu2[(
3255 nd + ii * sotbu2_dim2) * sotbu2_dim1] * fpntab[nuroo
3256 + 1 + (ii << 1) * fpntab_dim1];
3257 sosotb[nd * sosotb_dim2 * sosotb_dim1] -= bid1;
3266 /* ------------------------------ The End -------------------------------
3271 AdvApp2Var_SysBase::mgsomsg_("MMA2CD3", 7L);
3276 //=======================================================================
3277 //function : mma2cdi_
3279 //=======================================================================
3280 int AdvApp2Var_ApproxF2var::mma2cdi_( integer *ndimen,
3308 /* System generated locals */
3309 integer contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
3310 contr2_offset, contr3_dim1, contr3_dim2, contr3_offset,
3311 contr4_dim1, contr4_dim2, contr4_offset, sosotb_dim1, sosotb_dim2,
3312 sosotb_offset, diditb_dim1, diditb_dim2, diditb_offset,
3313 soditb_dim1, soditb_dim2, soditb_offset, disotb_dim1, disotb_dim2,
3316 /* Local variables */
3319 doublereal* wrkar = 0;
3321 integer ibb, ier = 0;
3322 integer isz1, isz2, isz3, isz4;
3323 intptr_t ipt1, ipt2, ipt3, ipt4;
3328 /* **********************************************************************
3333 /* Discretisation on the parameters of polynomes of interpolation */
3334 /* of constraints of order IORDRE. */
3338 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
3340 //* INPUT ARGUMENTS : */
3341 /* ------------------ */
3342 /* NDIMEN: Dimension of the space. */
3343 /* NBPNTU: Nb of INTERNAL parameters of discretisation by U. */
3344 /* This is also the nb of root of Legendre polynom where discretization is done. */
3345 /* UROOTL: Table of parameters of discretisation ON (-1,1) by U.
3347 /* NBPNTV: Nb of INTERNAL parameters of discretisation by V. */
3348 /* This is also the nb of root of Legendre polynom where discretization is done. */
3349 /* VROOTL: Table of parameters of discretisation ON (-1,1) by V.*/
3351 /* IORDRV: Order of constraint imposed at the extremities of iso-U */
3352 /* = 0, calculate the extremities of iso-U */
3353 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-U */
3354 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-U */
3355 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
3356 /* = 0, calculate the extremities of iso-V */
3357 /* = 1, calculate, additionally, the 1st derivative in the direction of iso-V */
3358 /* = 2, calculate, additionally, the 2nd derivative in the direction of iso-V */
3359 /* CONTR1: Contains, if IORDRU and IORDRV>=0, the values at the */
3360 /* extremities of F(U0,V0) and its derivatives. */
3361 /* CONTR2: Contains, if IORDRU and IORDRV>=0, the values at the */
3362 /* extremities of F(U1,V0) and its derivatives. */
3363 /* CONTR3: Contains, if IORDRU and IORDRV>=0, the values at the */
3364 /* extremities of F(U0,V1) and its derivatives. */
3365 /* CONTR4: Contains, if IORDRU and IORDRV>=0, the values at the */
3366 /* extremities of F(U1,V1) and its derivatives. */
3367 /* SOTBU1: Table of NBPNTU/2 sums of 2 index points */
3368 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3369 /* SOTBU2: Table of NBPNTV/2 sums of 2 index points */
3370 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3371 /* DITBU1: Table of NBPNTU/2 differences of 2 index points */
3372 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V0. */
3373 /* DITBU2: Table of NBPNTU/2 differences of 2 index points */
3374 /* NBPNTU-II+1 and II, for II = 1, NBPNTU/2 on iso-V1. */
3375 /* SOTBV1: Table of NBPNTV/2 sums of 2 index points */
3376 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
3377 /* SOTBV2: Table of NBPNTV/2 sums of 2 index points */
3378 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
3379 /* DITBV1: Table of NBPNTV/2 differences of 2 index points */
3380 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V0. */
3381 /* DITBV2: Table of NBPNTV/2 differences of 2 index points */
3382 /* NBPNTV-II+1 and II, for II = 1, NBPNTV/2 on iso-V1. */
3383 /* SOSOTB: Preinitialized table (input/output argument). */
3384 /* DISOTB: Preinitialized table (input/output argument). */
3385 /* SODITB: Preinitialized table (input/output argument). */
3386 /* DIDITB: Preinitialized table (input/output argument) */
3388 /* ARGUMENTS DE SORTIE : */
3389 /* ------------------- */
3390 /* SOSOTB: Table where the terms of constraints are added */
3391 /* C(ui,vj) + C(ui,-vj) + C(-ui,vj) + C(-ui,-vj) */
3392 /* with ui and vj positive roots of the Legendre polynom */
3393 /* of degree NBPNTU and NBPNTV respectively. */
3394 /* DISOTB: Table where the terms of constraints are added */
3395 /* C(ui,vj) + C(ui,-vj) - C(-ui,vj) - C(-ui,-vj) */
3396 /* with ui and vj positive roots of the polynom of Legendre */
3397 /* of degree NBPNTU and NBPNTV respectively. */
3398 /* SODITB: Table where the terms of constraints are added */
3399 /* C(ui,vj) - C(ui,-vj) + C(-ui,vj) - C(-ui,-vj) */
3400 /* with ui and vj positive roots of the polynom of Legendre */
3401 /* of degree NBPNTU and NBPNTV respectively. */
3402 /* DIDITB: Table where the terms of constraints are added */
3403 /* C(ui,vj) - C(ui,-vj) - C(-ui,vj) + C(-ui,-vj) */
3404 /* with ui and vj positive roots of the polynom of Legendre */
3405 /* of degree NBPNTU and NBPNTV respectively. */
3406 /* IERCOD: = 0, OK, */
3407 /* = 1, Value or IORDRV or IORDRU is out of allowed values. */
3408 /* =13, Pb of dynamic allocation. */
3410 /* COMMONS USED : */
3411 /* ---------------- */
3413 /* REFERENCES CALLED : */
3414 /* -------------------- */
3416 /* DESCRIPTION/NOTES/LIMITATIONS : */
3417 /* ------------------------------- */
3420 /* **********************************************************************
3423 /* The name of the routine */
3426 /* Parameter adjustments */
3428 diditb_dim1 = *nbpntu / 2 + 1;
3429 diditb_dim2 = *nbpntv / 2 + 1;
3430 diditb_offset = diditb_dim1 * diditb_dim2;
3431 diditb -= diditb_offset;
3432 disotb_dim1 = *nbpntu / 2;
3433 disotb_dim2 = *nbpntv / 2;
3434 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3435 disotb -= disotb_offset;
3436 soditb_dim1 = *nbpntu / 2;
3437 soditb_dim2 = *nbpntv / 2;
3438 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3439 soditb -= soditb_offset;
3440 sosotb_dim1 = *nbpntu / 2 + 1;
3441 sosotb_dim2 = *nbpntv / 2 + 1;
3442 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3443 sosotb -= sosotb_offset;
3445 contr4_dim1 = *ndimen;
3446 contr4_dim2 = *iordru + 2;
3447 contr4_offset = contr4_dim1 * (contr4_dim2 + 1) + 1;
3448 contr4 -= contr4_offset;
3449 contr3_dim1 = *ndimen;
3450 contr3_dim2 = *iordru + 2;
3451 contr3_offset = contr3_dim1 * (contr3_dim2 + 1) + 1;
3452 contr3 -= contr3_offset;
3453 contr2_dim1 = *ndimen;
3454 contr2_dim2 = *iordru + 2;
3455 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
3456 contr2 -= contr2_offset;
3457 contr1_dim1 = *ndimen;
3458 contr1_dim2 = *iordru + 2;
3459 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
3460 contr1 -= contr1_offset;
3469 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
3472 ibb = AdvApp2Var_SysBase::mnfndeb_();
3474 AdvApp2Var_SysBase::mgenmsg_("MMA2CDI", 7L);
3478 if (*iordru < -1 || *iordru > 2) {
3481 if (*iordrv < -1 || *iordrv > 2) {
3485 /* ------------------------- Set to zero --------------------------------
3488 ilong = (*nbpntu / 2 + 1) * (*nbpntv / 2 + 1) * *ndimen;
3489 AdvApp2Var_SysBase::mvriraz_(&ilong, &sosotb[sosotb_offset]);
3490 AdvApp2Var_SysBase::mvriraz_(&ilong, &diditb[diditb_offset]);
3491 ilong = *nbpntu / 2 * (*nbpntv / 2) * *ndimen;
3492 AdvApp2Var_SysBase::mvriraz_(&ilong, &soditb[soditb_offset]);
3493 AdvApp2Var_SysBase::mvriraz_(&ilong, &disotb[disotb_offset]);
3494 if (*iordru == -1 && *iordrv == -1) {
3500 isz1 = ((*iordru + 1) << 2) * (*iordru + 1);
3501 isz2 = ((*iordrv + 1) << 2) * (*iordrv + 1);
3502 isz3 = ((*iordru + 1) << 1) * *nbpntu;
3503 isz4 = ((*iordrv + 1) << 1) * *nbpntv;
3504 iszwr = isz1 + isz2 + isz3 + isz4;
3505 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3514 if (*iordru >= 0 && *iordru <= 2) {
3516 /* --- Return 2*(IORDRU+1) coeff of 2*(IORDRU+1) polynoms of Hermite
3519 AdvApp2Var_ApproxF2var::mma1her_(iordru, &wrkar[ipt1], iercod);
3524 /* ---- Subract discretizations of polynoms of constraints
3527 mma2cd3_(ndimen, nbpntu, &urootl[1], nbpntv, iordru, &sotbu1[1], &
3528 sotbu2[1], &ditbu1[1], &ditbu2[1], &wrkar[ipt3], &wrkar[ipt1],
3529 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3530 disotb_offset], &diditb[diditb_offset]);
3533 if (*iordrv >= 0 && *iordrv <= 2) {
3535 /* --- Return 2*(IORDRV+1) coeff of 2*(IORDRV+1) polynoms of Hermite
3538 AdvApp2Var_ApproxF2var::mma1her_(iordrv, &wrkar[ipt2], iercod);
3543 /* ---- Subtract discretisations of polynoms of constraint
3546 mma2cd2_(ndimen, nbpntu, nbpntv, &vrootl[1], iordrv, &sotbv1[1], &
3547 sotbv2[1], &ditbv1[1], &ditbv2[1], &wrkar[ipt4], &wrkar[ipt2],
3548 &sosotb[sosotb_offset], &soditb[soditb_offset], &disotb[
3549 disotb_offset], &diditb[diditb_offset]);
3552 /* --------------- Subtract constraints of corners ----------------
3555 if (*iordru >= 0 && *iordrv >= 0) {
3556 mma2cd1_(ndimen, nbpntu, &urootl[1], nbpntv, &vrootl[1], iordru,
3557 iordrv, &contr1[contr1_offset], &contr2[contr2_offset], &
3558 contr3[contr3_offset], &contr4[contr4_offset], &wrkar[ipt3], &
3559 wrkar[ipt4], &wrkar[ipt1], &wrkar[ipt2], &sosotb[
3560 sosotb_offset], &soditb[soditb_offset], &disotb[disotb_offset]
3561 , &diditb[diditb_offset]);
3565 /* ------------------------------ The End -------------------------------
3567 /* --> IORDRE is not within the autorised diapason. */
3571 /* --> PB of dynamic allocation. */
3578 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3583 AdvApp2Var_SysBase::maermsg_("MMA2CDI", iercod, 7L);
3585 AdvApp2Var_SysBase::mgsomsg_("MMA2CDI", 7L);
3590 //=======================================================================
3591 //function : mma2ce1_
3593 //=======================================================================
3594 int AdvApp2Var_ApproxF2var::mma2ce1_(integer *numdec,
3624 /* System generated locals */
3625 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3626 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3627 diditb_dim1, diditb_dim2, diditb_offset, patjac_dim1, patjac_dim2,
3630 /* Local variables */
3633 doublereal* wrkar = 0;
3636 integer isz1, isz2, isz3, isz4, isz5, isz6, isz7;
3637 intptr_t ipt1, ipt2, ipt3, ipt4, ipt5, ipt6, ipt7;
3641 /* **********************************************************************
3646 /* Calculation of coefficients of polynomial approximation of degree */
3647 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3648 /* discretization on roots of Legendre polynom of degree */
3649 /* NBPNTU by U and NBPNTV by V. */
3653 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&POLYNOME,&ERREUR */
3655 /* INPUT ARGUMENTS : */
3656 /* ------------------ */
3657 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3658 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3659 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3660 /* directions simultaneously (cutting by V is preferable). */
3661 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3662 /* directions simultaneously (cutting by U is preferable). */
3663 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3664 /* of cutting Vj). */
3665 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3666 /* of cutting Ui). */
3667 /* = 0, It is not POSSIBLE to cut anything */
3668 /* NDIMEN: Dimension of the space. */
3669 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3670 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3671 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3672 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3673 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3674 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3675 /* NDJACU: Max degree of the polynom of approximation by U. */
3676 /* The representation in the orthogonal base starts from degree */
3677 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3678 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3679 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3680 /* NDJACV: Max degree of the polynom of approximation by V. */
3681 /* The representation in the orthogonal base starts from degree */
3682 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3683 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3684 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3685 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3686 /* to the step of constraints C0, C1 or C2. */
3687 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3688 /* to the step of constraints C0, C1 or C2. */
3689 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3690 /* calculated the coefficients of integration by u */
3691 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3692 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3693 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3694 /* calculated the coefficients of integration by u */
3695 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3696 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3697 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3698 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3699 /* with ui and vj - positive roots of the Legendre polynom */
3700 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3701 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3702 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3703 /* SOSOTB(0,0) contains F(0,0). */
3704 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3705 /* with ui and vj positive roots of Legendre polynom */
3706 /* of degree NBPNTU and NBPNTV respectively. */
3707 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3708 /* with ui and vj positive roots of Legendre polynom */
3709 /* of degree NBPNTU and NBPNTV respectively. */
3710 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3711 /* with ui and vj positive roots of Legendre polynom */
3712 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3713 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3714 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3716 /* OUTPUT ARGUMENTS */
3717 /* --------------- */
3718 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
3719 /* of F(u,v) with eventually taking into account of */
3720 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
3721 /* This table contains other coeff if ITYDEC = 0. */
3722 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
3723 /* on each of sub-spaces SI ITYDEC = 0. */
3724 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
3725 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
3726 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
3727 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
3728 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
3729 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
3730 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
3731 /* = 3, it is NECESSARY to cut both by U AND by V. */
3732 /* IERCOD: Error code. */
3733 /* = 0, Everything is OK. */
3734 /* = -1, There is the best possible solution, but the */
3735 /* user tolerance is not satisfactory (3*only) */
3736 /* = 1, Incoherent entries. */
3738 /* COMMONS USED : */
3739 /* ---------------- */
3741 /* REFERENCES CALLED : */
3742 /* --------------------- */
3744 /* DESCRIPTION/NOTES/LIMITATIONS : */
3745 /* ------------------------------- */
3748 /* **********************************************************************
3750 /* Name of the routine */
3753 /* --------------------------- Initialisations --------------------------
3756 /* Parameter adjustments */
3761 patjac_dim1 = *ndjacu + 1;
3762 patjac_dim2 = *ndjacv + 1;
3763 patjac_offset = patjac_dim1 * patjac_dim2;
3764 patjac -= patjac_offset;
3765 diditb_dim1 = *nbpntu / 2 + 1;
3766 diditb_dim2 = *nbpntv / 2 + 1;
3767 diditb_offset = diditb_dim1 * diditb_dim2;
3768 diditb -= diditb_offset;
3769 soditb_dim1 = *nbpntu / 2;
3770 soditb_dim2 = *nbpntv / 2;
3771 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
3772 soditb -= soditb_offset;
3773 disotb_dim1 = *nbpntu / 2;
3774 disotb_dim2 = *nbpntv / 2;
3775 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
3776 disotb -= disotb_offset;
3777 sosotb_dim1 = *nbpntu / 2 + 1;
3778 sosotb_dim2 = *nbpntv / 2 + 1;
3779 sosotb_offset = sosotb_dim1 * sosotb_dim2;
3780 sosotb -= sosotb_offset;
3783 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
3785 AdvApp2Var_SysBase::mgenmsg_("MMA2CE1", 7L);
3790 isz1 = (*nbpntu / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1);
3791 isz2 = (*nbpntv / 2 + 1) * (*ndjacv - ((*iordrv + 1) << 1) + 1);
3792 isz3 = (*nbpntv / 2 + 1) * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3793 isz4 = *nbpntv / 2 * (*ndjacu - ((*iordru + 1) << 1) + 1) * *ndimen;
3794 isz5 = *ndjacu + 1 - ((*iordru + 1) << 1);
3795 isz6 = *ndjacv + 1 - ((*iordrv + 1) << 1);
3796 isz7 = *ndimen << 2;
3797 iszwr = isz1 + isz2 + isz3 + isz4 + isz5 + isz6 + isz7;
3798 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
3799 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3811 /* ----------------- Return Gauss coefficients of integration ----------------
3814 AdvApp2Var_ApproxF2var::mmapptt_(ndjacu, nbpntu, iordru, &wrkar[ipt1], iercod);
3818 AdvApp2Var_ApproxF2var::mmapptt_(ndjacv, nbpntv, iordrv, &wrkar[ipt2], iercod);
3823 /* ------------------- Return max polynoms of Jacobi ------------
3826 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacu, iordru, &wrkar[ipt5]);
3827 AdvApp2Var_ApproxF2var::mma2jmx_(ndjacv, iordrv, &wrkar[ipt6]);
3829 /* ------ Calculate the coefficients and their contribution to the error ----
3832 mma2ce2_(numdec, ndimen, nbsesp, &ndimse[1], ndminu, ndminv, ndguli,
3833 ndgvli, ndjacu, ndjacv, iordru, iordrv, nbpntu, nbpntv, &epsapr[1]
3834 , &sosotb[sosotb_offset], &disotb[disotb_offset], &soditb[
3835 soditb_offset], &diditb[diditb_offset], &wrkar[ipt1], &wrkar[ipt2]
3836 , &wrkar[ipt5], &wrkar[ipt6], &wrkar[ipt7], &wrkar[ipt3], &wrkar[
3837 ipt4], &patjac[patjac_offset], &errmax[1], &errmoy[1], ndegpu,
3838 ndegpv, itydec, iercod);
3844 /* ------------------------------ The end -------------------------------
3853 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &iszwr, wrkar, &iofwr, &ier);
3858 AdvApp2Var_SysBase::maermsg_("MMA2CE1", iercod, 7L);
3860 AdvApp2Var_SysBase::mgsomsg_("MMA2CE1", 7L);
3865 //=======================================================================
3866 //function : mma2ce2_
3868 //=======================================================================
3869 int mma2ce2_(integer *numdec,
3904 /* System generated locals */
3905 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
3906 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
3907 diditb_dim1, diditb_dim2, diditb_offset, gssutb_dim1, gssvtb_dim1,
3908 chpair_dim1, chpair_dim2, chpair_offset, chimpr_dim1,
3909 chimpr_dim2, chimpr_offset, patjac_dim1, patjac_dim2,
3910 patjac_offset, vecerr_dim1, vecerr_offset, i__1, i__2, i__3, i__4;
3912 /* Local variables */
3914 integer idim, igsu, minu, minv, maxu, maxv, igsv;
3916 integer i2rdu, i2rdv, ndses, nd, ii, jj, kk, nu, nv;
3920 /* **********************************************************************
3924 /* Calculation of coefficients of polynomial approximation of degree */
3925 /* (NDJACU,NDJACV) of a function F(u,v), starting from its */
3926 /* discretization on roots of Legendre polynom of degree */
3927 /* NBPNTU by U and NBPNTV by V. */
3931 /* TOUS,AB_SPECIFI::FONCTION&,APPROXIMATION,&COEFFICIENT,&POLYNOME */
3933 /* INPUT ARGUMENTS : */
3934 /* ------------------ */
3935 /* NUMDEC: Indicates if it is POSSIBLE to cut function F(u,v). */
3936 /* = 5, It is POSSIBLE to cut by U or by V or in both directions simultaneously. */
3937 /* = 4, It is POSSIBLE to cut by U or by V BUT NOT in both */
3938 /* directions simultaneously (cutting by V is preferable). */
3939 /* = 3, It is POSSIBLE to cut by U or by V BUT NOT in both */
3940 /* directions simultaneously (cutting by U is preferable). */
3941 /* = 2, It is POSSIBLE to cut only by V (i.e. insert parameter */
3942 /* of cutting Vj). */
3943 /* = 1, It is POSSIBLE to cut only by U (i.e. insert parameter */
3944 /* of cutting Ui). */
3945 /* = 0, It is not POSSIBLE to cut anything */
3946 /* NDIMEN: Total dimension of the space. */
3947 /* NBSESP: Nb of independent sub-spaces on which the errors are calculated. */
3948 /* NDIMSE: Table of dimensions of each of sub-spaces. */
3949 /* NDMINU: Minimum degree by U to be preserved for the approximation. */
3950 /* NDMINV: Minimum degree by V to be preserved for the approximation. */
3951 /* NDGULI: Limit of nb of coefficients by U of the solution. */
3952 /* NDGVLI: Limit of nb of coefficients by V of the solution. */
3953 /* NDJACU: Max degree of the polynom of approximation by U. */
3954 /* The representation in the orthogonal base starts from degree */
3955 /* 0 to degree NDJACU-2*(IORDRU+1). The polynomial base is the base of */
3956 /* Jacobi of order -1 (Legendre), 0, 1 or 2. */
3957 /* It is required that 2*IORDRU+1 <= NDMINU <= NDGULI < NDJACU */
3958 /* NDJACV: Max degree of the polynom of approximation by V. */
3959 /* The representation in the orthogonal base starts from degree */
3960 /* 0 to degree NDJACV-2*(IORDRV+1). The polynomial base is */
3961 /* the base of Jacobi of order -1 (Legendre), 0, 1 or 2 */
3962 /* It is required that 2*IORDRV+1 <= NDMINV <= NDGVLI < NDJACV */
3963 /* IORDRU: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3964 /* to the step of constraints C0, C1 or C2. */
3965 /* IORDRV: Order of the Jacobi base (-1,0,1 or 2) by U. Corresponds */
3966 /* to the step of constraints C0, C1 or C2. */
3967 /* NBPNTU: Degree of Legendre polynom on the roots which of are */
3968 /* calculated the coefficients of integration by u */
3969 /* by Gauss method. It is required that NBPNTU = 30, 40, */
3970 /* 50 or 61 and NDJACU-2*(IORDRU+1) < NBPNTU. */
3971 /* NBPNTV: Degree of Legendre polynom on the roots which of are */
3972 /* calculated the coefficients of integration by u */
3973 /* by Gauss method. It is required that NBPNTV = 30, 40, */
3974 /* 50 or 61 and NDJACV-2*(IORDRV+1) < NBPNTV. */
3975 /* EPSAPR: Table of NBSESP tolerances imposed on each sub-spaces. */
3976 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
3977 /* with ui and vj - positive roots of the Legendre polynom */
3978 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3979 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
3980 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
3981 /* SOSOTB(0,0) contains F(0,0). */
3982 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
3983 /* with ui and vj positive roots of Legendre polynom */
3984 /* of degree NBPNTU and NBPNTV respectively. */
3985 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
3986 /* with ui and vj positive roots of Legendre polynom */
3987 /* of degree NBPNTU and NBPNTV respectively. */
3988 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
3989 /* with ui and vj positive roots of Legendre polynom */
3990 /* of degree NBPNTU and NBPNTV respectively. Additionally, */
3991 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
3992 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
3993 /* GSSUTB: Table of coefficients of integration by Gauss method */
3994 /* by U: i varies from 0 to NBPNTU/2 and k varies from 0 to */
3995 /* NDJACU-2*(IORDRU+1). */
3996 /* GSSVTB: Table of coefficients of integration by Gauss method */
3997 /* by V: i varies from 0 to NBPNTV/2 and k varies from 0 to */
3998 /* NDJACV-2*(IORDRV+1). */
3999 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
4000 /* from degree 0 to degree NDJACU - 2*(IORDRU+1) */
4001 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
4002 /* from degree 0 to degree NDJACV - 2*(IORDRV+1) */
4004 /* OUTPUT ARGUMENTS : */
4005 /* ------------------- */
4006 /* VECERR: Auxiliary table. */
4007 /* CHPAIR: Auxiliary table of terms connected to degree NDJACU by U */
4008 /* to calculate the coeff. of approximation of EVEN degree by V. */
4009 /* CHIMPR: Auxiliary table of terms connected to degree NDJACU by U */
4010 /* to calculate the coeff. of approximation of UNEVEN degree by V. */
4011 /* PATJAC: Table of coefficients of polynom P(u,v) of approximation */
4012 /* of F(u,v) with eventually taking into account of */
4013 /* constraints. P(u,v) is of degree (NDJACU,NDJACV). */
4014 /* This table contains other coeff if ITYDEC = 0. */
4015 /* ERRMAX: For 1<=i<=NBSESP, ERRMAX(i) contains max errors */
4016 /* on each of sub-spaces SI ITYDEC = 0. */
4017 /* ERRMOY: Contains average errors for each of NBSESP sub-spaces SI ITYDEC = 0. */
4018 /* NDEGPU: Degree by U for square PATJAC. Valable if ITYDEC=0. */
4019 /* NDEGPV: Degree by V for square PATJAC. Valable if ITYDEC=0. */
4020 /* ITYDEC: Shows if it is NECESSARY to cut again function F(u,v). */
4021 /* = 0, it is not NECESSARY to cut anything, PATJAC is OK. */
4022 /* = 1, it is NECESSARY to cut only by U (i.e. insert parameter of cutting Ui). */
4023 /* = 2, it is NECESSARY to cut only by V (i.e. insert parameter of cutting Vj). */
4024 /* = 3, it is NECESSARY to cut both by U AND by V. */
4025 /* IERCOD: Error code. */
4026 /* = 0, Everything is OK. */
4027 /* = -1, There is the best possible solution, but the */
4028 /* user tolerance is not satisfactory (3*only) */
4029 /* = 1, Incoherent entries. */
4031 /* COMMONS USED : */
4032 /* ---------------- */
4034 /* REFERENCES CALLED : */
4035 /* --------------------- */
4037 /* DESCRIPTION/NOTES/LIMITATIONS : */
4039 /* **********************************************************************
4041 /* Name of the routine */
4044 /* --------------------------- Initialisations --------------------------
4047 /* Parameter adjustments */
4048 vecerr_dim1 = *ndimen;
4049 vecerr_offset = vecerr_dim1 + 1;
4050 vecerr -= vecerr_offset;
4055 patjac_dim1 = *ndjacu + 1;
4056 patjac_dim2 = *ndjacv + 1;
4057 patjac_offset = patjac_dim1 * patjac_dim2;
4058 patjac -= patjac_offset;
4059 gssutb_dim1 = *nbpntu / 2 + 1;
4060 chimpr_dim1 = *nbpntv / 2;
4061 chimpr_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4062 chimpr_offset = chimpr_dim1 * chimpr_dim2 + 1;
4063 chimpr -= chimpr_offset;
4064 chpair_dim1 = *nbpntv / 2 + 1;
4065 chpair_dim2 = *ndjacu - ((*iordru + 1) << 1) + 1;
4066 chpair_offset = chpair_dim1 * chpair_dim2;
4067 chpair -= chpair_offset;
4068 gssvtb_dim1 = *nbpntv / 2 + 1;
4069 diditb_dim1 = *nbpntu / 2 + 1;
4070 diditb_dim2 = *nbpntv / 2 + 1;
4071 diditb_offset = diditb_dim1 * diditb_dim2;
4072 diditb -= diditb_offset;
4073 soditb_dim1 = *nbpntu / 2;
4074 soditb_dim2 = *nbpntv / 2;
4075 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
4076 soditb -= soditb_offset;
4077 disotb_dim1 = *nbpntu / 2;
4078 disotb_dim2 = *nbpntv / 2;
4079 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
4080 disotb -= disotb_offset;
4081 sosotb_dim1 = *nbpntu / 2 + 1;
4082 sosotb_dim2 = *nbpntv / 2 + 1;
4083 sosotb_offset = sosotb_dim1 * sosotb_dim2;
4084 sosotb -= sosotb_offset;
4087 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4089 AdvApp2Var_SysBase::mgenmsg_("MMA2CE2", 7L);
4091 /* --> A priori everything is OK */
4093 /* --> test of inputs */
4094 if (*numdec < 0 || *numdec > 5) {
4097 if ((*iordru << 1) + 1 > *ndminu) {
4100 if (*ndminu > *ndguli) {
4103 if (*ndguli >= *ndjacu) {
4106 if ((*iordrv << 1) + 1 > *ndminv) {
4109 if (*ndminv > *ndgvli) {
4112 if (*ndgvli >= *ndjacv) {
4115 /* --> A priori, no cuts to be done */
4117 /* --> Min. degrees to return: NDMINU,NDMINV */
4120 /* --> For the moment, max errors are null */
4121 AdvApp2Var_SysBase::mvriraz_(nbsesp, &errmax[1]);
4123 AdvApp2Var_SysBase::mvriraz_(&nd, &vecerr[vecerr_offset]);
4124 /* --> and the square, too. */
4125 nd = (*ndjacu + 1) * (*ndjacv + 1) * *ndimen;
4126 AdvApp2Var_SysBase::mvriraz_(&nd, &patjac[patjac_offset]);
4128 i2rdu = (*iordru + 1) << 1;
4129 i2rdv = (*iordrv + 1) << 1;
4131 /* **********************************************************************
4133 /* -------------------- HERE IT IS POSSIBLE TO CUT ----------------------
4135 /* **********************************************************************
4138 if (*numdec > 0 && *numdec <= 5) {
4140 /* ******************************************************************
4142 /* ---------------------- Calculate coeff of zone 4 -------------
4156 /* ---------------- Calculate the terms connected to degree by U ---------
4160 for (nd = 1; nd <= i__1; ++nd) {
4162 for (kk = minu; kk <= i__2; ++kk) {
4164 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4165 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4166 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4167 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4168 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4169 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4170 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4176 /* ------------------- Calculate the coefficients of PATJAC ------------
4179 igsu = minu - i2rdu;
4181 for (jj = minv; jj <= i__1; ++jj) {
4184 for (nd = 1; nd <= i__2; ++nd) {
4185 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4186 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4187 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4188 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4189 patjac_dim2) * patjac_dim1]);
4193 /* ----- Contribution of calculated terms to the approximation error */
4194 /* for terms (I,J) with MINU <= I <= MAXU, J fixe. */
4198 for (nd = 1; nd <= i__2; ++nd) {
4200 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &jj, &jj,
4201 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4202 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4203 &vecerr[nd + (vecerr_dim1 << 2)]);
4204 if (vecerr[nd + (vecerr_dim1 << 2)] > epsapr[nd]) {
4213 /* ******************************************************************
4215 /* ---------------------- Calculate the coeff of zone 2 -------------
4218 minu = (*iordru + 1) << 1;
4223 /* --> If zone 2 is empty, pass to zone 3. */
4224 /* VECERR(ND,2) was already set to zero. */
4229 /* ---------------- Calculate the terms connected to degree by U ------------
4233 for (nd = 1; nd <= i__1; ++nd) {
4235 for (kk = minu; kk <= i__2; ++kk) {
4237 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4238 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4239 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4240 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4241 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4242 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4243 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4249 /* ------------------- Calculate the coefficients of PATJAC ------------
4252 igsu = minu - i2rdu;
4254 for (jj = minv; jj <= i__1; ++jj) {
4257 for (nd = 1; nd <= i__2; ++nd) {
4258 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4259 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4260 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4261 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4262 patjac_dim2) * patjac_dim1]);
4268 /* -----Contribution of calculated terms to the approximation error */
4269 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4273 for (nd = 1; nd <= i__1; ++nd) {
4275 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4276 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4277 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4278 vecerr[nd + (vecerr_dim1 << 1)]);
4283 /* ******************************************************************
4285 /* ---------------------- Calculation of coeff of zone 3 -------------
4291 minv = (*iordrv + 1) << 1;
4294 /* -> If zone 3 is empty, pass to the test of cutting. */
4295 /* VECERR(ND,3) was already set to zero */
4300 /* ----------- The terms connected to the degree by U are already calculated -----
4302 /* ------------------- Calculation of coefficients of PATJAC ------------
4305 igsu = minu - i2rdu;
4307 for (jj = minv; jj <= i__1; ++jj) {
4310 for (nd = 1; nd <= i__2; ++nd) {
4311 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4312 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4313 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4314 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4315 patjac_dim2) * patjac_dim1]);
4321 /* ----- Contribution of calculated terms to the approximation error */
4322 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV. */
4326 for (nd = 1; nd <= i__1; ++nd) {
4328 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4329 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4330 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4331 vecerr[nd + vecerr_dim1 * 3]);
4336 /* ******************************************************************
4338 /* --------------------------- Tests of cutting ---------------------
4343 for (nd = 1; nd <= i__1; ++nd) {
4344 vaux[0] = vecerr[nd + (vecerr_dim1 << 1)];
4345 vaux[1] = vecerr[nd + (vecerr_dim1 << 2)];
4346 vaux[2] = vecerr[nd + vecerr_dim1 * 3];
4348 errmax[nd] = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4349 if (errmax[nd] > epsapr[nd]) {
4351 zv = AdvApp2Var_MathBase::mzsnorm_(&ii, vaux);
4352 zu = AdvApp2Var_MathBase::mzsnorm_(&ii, &vaux[1]);
4353 if (zu > epsapr[nd] && zv > epsapr[nd]) {
4365 /* ******************************************************************
4367 /* --- OK, the square is valid, the coeff of zone 1 are calculated
4370 minu = (*iordru + 1) << 1;
4372 minv = (*iordrv + 1) << 1;
4375 /* --> If zone 1 is empty, pass to the calculation of Max and Average error. */
4376 if (minu > maxu || minv > maxv) {
4380 /* ----------- The terms connected to degree by U are already calculated -----
4382 /* ------------------- Calculate the coefficients of PATJAC ------------
4385 igsu = minu - i2rdu;
4387 for (jj = minv; jj <= i__1; ++jj) {
4390 for (nd = 1; nd <= i__2; ++nd) {
4391 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4392 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4393 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4394 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4395 patjac_dim2) * patjac_dim1]);
4401 /* --------------- Now the degree is maximally lowered --------
4406 i__1 = 1, i__2 = (*iordru << 1) + 1, i__1 = advapp_max(i__1,i__2);
4407 minu = advapp_max(i__1,*ndminu);
4410 i__1 = 1, i__2 = (*iordrv << 1) + 1, i__1 = advapp_max(i__1,i__2);
4411 minv = advapp_max(i__1,*ndminv);
4415 for (nd = 1; nd <= i__1; ++nd) {
4417 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4418 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4419 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4420 patjac_dim2 * patjac_dim1], &epsapr[nd], &vecerr[
4421 vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4429 /* --> Calculate the average error. */
4430 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4431 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4434 /* --> Set to 0.D0 the rejected coeffs. */
4435 i__2 = idim + ndses - 1;
4436 for (ii = idim; ii <= i__2; ++ii) {
4438 for (jj = nv1; jj <= i__3; ++jj) {
4440 for (kk = nu1; kk <= i__4; ++kk) {
4441 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4450 /* --> Return the nb of coeffs of approximation. */
4451 *ndegpu = advapp_max(*ndegpu,nu);
4452 *ndegpv = advapp_max(*ndegpv,nv);
4457 /* ******************************************************************
4459 /* -------------------- IT IS NOT POSSIBLE TO CUT -------------------
4461 /* ******************************************************************
4465 minu = (*iordru + 1) << 1;
4467 minv = (*iordrv + 1) << 1;
4470 /* ---------------- Calculate the terms connected to the degree by U ------------
4474 for (nd = 1; nd <= i__1; ++nd) {
4476 for (kk = minu; kk <= i__2; ++kk) {
4478 mma2cfu_(&kk, nbpntu, nbpntv, &sosotb[nd * sosotb_dim2 *
4479 sosotb_dim1], &disotb[(nd * disotb_dim2 + 1) *
4480 disotb_dim1 + 1], &soditb[(nd * soditb_dim2 + 1) *
4481 soditb_dim1 + 1], &diditb[nd * diditb_dim2 *
4482 diditb_dim1], &gssutb[igsu * gssutb_dim1], &chpair[(
4483 igsu + nd * chpair_dim2) * chpair_dim1], &chimpr[(
4484 igsu + nd * chimpr_dim2) * chimpr_dim1 + 1]);
4488 /* ---------------------- Calculate all coefficients -------
4491 igsu = minu - i2rdu;
4493 for (jj = minv; jj <= i__2; ++jj) {
4495 mma2cfv_(&jj, &minu, &maxu, nbpntv, &gssvtb[igsv *
4496 gssvtb_dim1], &chpair[(igsu + nd * chpair_dim2) *
4497 chpair_dim1], &chimpr[(igsu + nd * chimpr_dim2) *
4498 chimpr_dim1 + 1], &patjac[minu + (jj + nd *
4499 patjac_dim2) * patjac_dim1]);
4505 /* ----- Contribution of calculated terms to the approximation error */
4506 /* for terms (I,J) with MINU <= I <= MAXU, MINV <= J <= MAXV */
4510 for (nd = 1; nd <= i__1; ++nd) {
4512 minu = (*iordru + 1) << 1;
4516 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4517 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4518 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1], &
4522 minv = (*iordrv + 1) << 1;
4525 mma2er1_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &maxv,
4526 iordru, iordrv, xmaxju, xmaxjv, &patjac[idim *
4527 patjac_dim2 * patjac_dim1], &vecerr[vecerr_dim1 + 1],
4531 /* ---------------------------- IF ERRMAX > EPSAPR, stop --------
4534 if (errmax[nd] > epsapr[nd]) {
4539 /* ------------- Otherwise, try to remove again the coeff
4544 i__2 = 1, i__3 = (*iordru << 1) + 1, i__2 = advapp_max(i__2,i__3);
4545 minu = advapp_max(i__2,*ndminu);
4548 i__2 = 1, i__3 = (*iordrv << 1) + 1, i__2 = advapp_max(i__2,i__3);
4549 minv = advapp_max(i__2,*ndminv);
4551 if (maxu >= (*iordru + 1) << 1 && maxv >= (*iordrv + 1) << 1) {
4552 mma2er2_(ndjacu, ndjacv, &ndses, &minu, &maxu, &minv, &
4553 maxv, iordru, iordrv, xmaxju, xmaxjv, &patjac[
4554 idim * patjac_dim2 * patjac_dim1], &epsapr[nd], &
4555 vecerr[vecerr_dim1 + 1], &errmax[nd], &nu, &nv);
4562 /* --------------------- Calculate the average error -------------
4567 mma2moy_(ndjacu, ndjacv, &ndses, &nu1, ndjacu, &nv1, ndjacv,
4568 iordru, iordrv, &patjac[idim * patjac_dim2 * patjac_dim1],
4571 /* --------------------- Set to 0.D0 the rejected coeffs ----------
4574 i__2 = idim + ndses - 1;
4575 for (ii = idim; ii <= i__2; ++ii) {
4577 for (jj = nv1; jj <= i__3; ++jj) {
4579 for (kk = nu1; kk <= i__4; ++kk) {
4580 patjac[kk + (jj + ii * patjac_dim2) * patjac_dim1] =
4589 /* --------------- Return the nb of coeff of approximation ---
4592 *ndegpu = advapp_max(*ndegpu,nu);
4593 *ndegpv = advapp_max(*ndegpv,nv);
4601 /* ------------------------------ The end -------------------------------
4603 /* --> Error in inputs */
4608 /* --------- Management of cuts, it is required 0 < NUMDEC <= 5 -------
4611 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4613 if (*numdec <= 0 || *numdec > 5) {
4622 /* --> Here it is possible and necessary to cut, choose by U if it is possible */
4624 if (*numdec <= 0 || *numdec > 5) {
4633 /* --> Here it is possible and necessary to cut, choose by 4 if it is possible */
4635 if (*numdec <= 0 || *numdec > 5) {
4640 } else if (*numdec == 2 || *numdec == 4) {
4642 } else if (*numdec == 1 || *numdec == 3) {
4650 AdvApp2Var_SysBase::maermsg_("MMA2CE2", iercod, 7L);
4652 AdvApp2Var_SysBase::mgsomsg_("MMA2CE2", 7L);
4657 //=======================================================================
4658 //function : mma2cfu_
4660 //=======================================================================
4661 int mma2cfu_(integer *ndujac,
4673 /* System generated locals */
4674 integer sosotb_dim1, disotb_dim1, disotb_offset, soditb_dim1,
4675 soditb_offset, diditb_dim1, i__1, i__2;
4677 /* Local variables */
4679 integer nptu2, nptv2, ii, jj;
4680 doublereal bid0, bid1, bid2;
4682 /* **********************************************************************
4687 /* Calculate the terms connected to degree NDUJAC by U of the polynomial approximation */
4688 /* of function F(u,v), starting from its discretisation */
4689 /* on the roots of Legendre polynom of degree */
4690 /* NBPNTU by U and NBPNTV by V. */
4694 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4696 /* INPUT ARGUMENTSE : */
4697 /* ------------------ */
4698 /* NDUJAC: Fixed degree by U for which the terms */
4699 /* allowing to obtain the Legendre or Jacobi coeff*/
4700 /* of even or uneven degree by V are calculated. */
4701 /* NBPNTU: Degree of Legendre polynom on the roots which of */
4702 /* the coefficients of integration by U are calculated */
4703 /* by Gauss method. It is required that NBPNTU = 30, 40, 50 or 61. */
4704 /* NBPNTV: Degree of Legendre polynom on the roots which of */
4705 /* the coefficients of integration by V are calculated */
4706 /* by Gauss method. It is required that NBPNTV = 30, 40, 50 or 61. */
4707 /* SOSOTB: Table of F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
4708 /* with ui and vj positive roots of Legendre polynom */
4709 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4710 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
4711 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0) and */
4712 /* SOSOTB(0,0) contains F(0,0). */
4713 /* DISOTB: Table of F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
4714 /* with ui and vj positive roots of Legendre polynom */
4715 /* of degree NBPNTU and NBPNTV respectively. */
4716 /* SODITB: Table of F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
4717 /* with ui and vj positive roots of Legendre polynom */
4718 /* of degree NBPNTU and NBPNTV respectively. */
4719 /* DIDITB: Table of F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
4720 /* avec ui and vj positive roots of Legendre polynom */
4721 /* of degree NBPNTU and NBPNTV respectively. Moreover, */
4722 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
4723 /* and table DIDITB(i,0) contains F(ui,0) - F(-ui,0). */
4724 /* GSSUTB: Table of coefficients of integration by Gauss method */
4725 /* Gauss by U for fixed NDUJAC : i varies from 0 to NBPNTU/2. */
4727 /* OUTPUT ARGUMENTS : */
4728 /* ------------------- */
4729 /* CHPAIR: Table of terms connected to degree NDUJAC by U to calculate the */
4730 /* coeff. of the approximation of EVEN degree by V. */
4731 /* CHIMPR: Table of terms connected to degree NDUJAC by U to calculate */
4732 /* the coeff. of approximation of UNEVEN degree by V. */
4734 /* COMMONS USED : */
4735 /* ---------------- */
4737 /* REFERENCES CALLED : */
4738 /* ----------------------- */
4740 /* DESCRIPTION/NOTES/LIMITATIONS : */
4741 /* ----------------------------------- */
4745 /* **********************************************************************
4747 /* Name of the routine */
4750 /* --------------------------- Initialisations --------------------------
4753 /* Parameter adjustments */
4755 diditb_dim1 = *nbpntu / 2 + 1;
4756 soditb_dim1 = *nbpntu / 2;
4757 soditb_offset = soditb_dim1 + 1;
4758 soditb -= soditb_offset;
4759 disotb_dim1 = *nbpntu / 2;
4760 disotb_offset = disotb_dim1 + 1;
4761 disotb -= disotb_offset;
4762 sosotb_dim1 = *nbpntu / 2 + 1;
4765 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4767 AdvApp2Var_SysBase::mgenmsg_("MMA2CFU", 7L);
4770 nptu2 = *nbpntu / 2;
4771 nptv2 = *nbpntv / 2;
4773 /* **********************************************************************
4775 /* CALCULATE COEFFICIENTS BY U */
4777 /* ----------------- Calculate coefficients of even degree --------------
4780 if (*ndujac % 2 == 0) {
4782 for (jj = 1; jj <= i__1; ++jj) {
4786 for (ii = 1; ii <= i__2; ++ii) {
4788 bid1 += sosotb[ii + jj * sosotb_dim1] * bid0;
4789 bid2 += soditb[ii + jj * soditb_dim1] * bid0;
4797 /* --------------- Calculate coefficients of uneven degree ----------
4802 for (jj = 1; jj <= i__1; ++jj) {
4806 for (ii = 1; ii <= i__2; ++ii) {
4808 bid1 += disotb[ii + jj * disotb_dim1] * bid0;
4809 bid2 += diditb[ii + jj * diditb_dim1] * bid0;
4818 /* ------- Add terms connected to the supplementary root (0.D0) ------ */
4819 /* ----------- of Legendre polynom of uneven degree NBPNTU -----------
4821 /* --> Only even NDUJAC terms are modified as GSSUTB(0) = 0 */
4822 /* when NDUJAC is uneven. */
4824 if (*nbpntu % 2 != 0 && *ndujac % 2 == 0) {
4827 for (jj = 1; jj <= i__1; ++jj) {
4828 chpair[jj] += sosotb[jj * sosotb_dim1] * bid0;
4829 chimpr[jj] += diditb[jj * diditb_dim1] * bid0;
4834 /* ------ Calculate the terms connected to supplementary roots (0.D0) ------
4836 /* ----------- of Legendre polynom of uneven degree NBPNTV -----------
4839 if (*nbpntv % 2 != 0) {
4840 /* --> Only CHPAIR terms are calculated as GSSVTB(0,IH-IDEBV)=0
4842 /* when IH is uneven (see MMA2CFV). */
4844 if (*ndujac % 2 == 0) {
4847 for (ii = 1; ii <= i__1; ++ii) {
4848 bid1 += sosotb[ii] * gssutb[ii];
4855 for (ii = 1; ii <= i__1; ++ii) {
4856 bid1 += diditb[ii] * gssutb[ii];
4861 if (*nbpntu % 2 != 0) {
4862 chpair[0] += sosotb[0] * gssutb[0];
4866 /* ------------------------------ The end -------------------------------
4870 AdvApp2Var_SysBase::mgsomsg_("MMA2CFU", 7L);
4875 //=======================================================================
4876 //function : mma2cfv_
4878 //=======================================================================
4879 int mma2cfv_(integer *ndvjac,
4889 /* System generated locals */
4890 integer chpair_dim1, chpair_offset, chimpr_dim1, chimpr_offset,
4891 patjac_offset, i__1, i__2;
4893 /* Local variables */
4895 integer nptv2, ii, jj;
4898 /* **********************************************************************
4903 /* Calculate the coefficients of polynomial approximation of F(u,v) */
4904 /* of degree NDVJAC by V and of degree by U varying from MINDGU to MAXDGU.
4909 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
4911 /* INPUT ARGUMENTS : */
4912 /* ------------------ */
4914 /* NDVJAC: Degree of the polynom of approximation by V. */
4915 /* The representation in the orthogonal base starts from degre 0. */
4916 /* The polynomial base is the base of Jacobi of order -1 */
4917 /* (Legendre), 0, 1 or 2 */
4918 /* MINDGU: Degree minimum by U of coeff. to calculate. */
4919 /* MAXDGU: Degree maximum by U of coeff. to calculate. */
4920 /* NBPNTV: Degree of the Legendre polynom on the roots which of */
4921 /* the coefficients of integration by V are calculated */
4922 /* by Gauss method. It is reqired that NBPNTV = 30, 40, 50 or 61 and NDVJAC < NBPNTV. */
4923 /* GSSVTB: Table of coefficients of integration by Gauss method */
4924 /* by V for NDVJAC fixed: j varies from 0 to NBPNTV/2. */
4925 /* CHPAIR: Table of terms connected to degrees from MINDGU to MAXDGU by U to */
4926 /* calculate the coeff. of approximation of EVEN degree NDVJAC by V. */
4927 /* CHIMPR: Table of terms connected to degrees from MINDGU to MAXDGU by U to */
4928 /* calculate the coeff. of approximation of UNEVEN degree NDVJAC by V. */
4930 /* OUTPUT ARGUMENTS : */
4931 /* ------------------- */
4932 /* PATJAC: Table of coefficients by U of the polynom of approximation */
4933 /* P(u,v) of degree MINDGU to MAXDGU by U and NDVJAC by V. */
4935 /* COMMONS USED : */
4936 /* -------------- */
4938 /* REFERENCES CALLED : */
4939 /* --------------------- */
4941 /* DESCRIPTION/NOTES/LIMITATIONS : */
4942 /* ------------------------------- */
4944 /* **********************************************************************
4946 /* Name of the routine */
4949 /* --------------------------- Initialisations --------------------------
4952 /* Parameter adjustments */
4953 patjac_offset = *mindgu;
4954 patjac -= patjac_offset;
4955 chimpr_dim1 = *nbpntv / 2;
4956 chimpr_offset = chimpr_dim1 * *mindgu + 1;
4957 chimpr -= chimpr_offset;
4958 chpair_dim1 = *nbpntv / 2 + 1;
4959 chpair_offset = chpair_dim1 * *mindgu;
4960 chpair -= chpair_offset;
4963 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
4965 AdvApp2Var_SysBase::mgenmsg_("MMA2CFV", 7L);
4967 nptv2 = *nbpntv / 2;
4969 /* --------- Calculate the coefficients for even degree NDVJAC ----------
4972 if (*ndvjac % 2 == 0) {
4974 for (ii = *mindgu; ii <= i__1; ++ii) {
4977 for (jj = 1; jj <= i__2; ++jj) {
4978 bid1 += chpair[jj + ii * chpair_dim1] * gssvtb[jj];
4985 /* -------- Calculate the coefficients for uneven degree NDVJAC -----
4990 for (ii = *mindgu; ii <= i__1; ++ii) {
4993 for (jj = 1; jj <= i__2; ++jj) {
4994 bid1 += chimpr[jj + ii * chimpr_dim1] * gssvtb[jj];
5002 /* ------- Add terms connected to the supplementary root (0.D0) ----- */
5003 /* --------of the Legendre polynom of uneven degree NBPNTV --------- */
5005 if (*nbpntv % 2 != 0 && *ndvjac % 2 == 0) {
5008 for (ii = *mindgu; ii <= i__1; ++ii) {
5009 patjac[ii] += bid1 * chpair[ii * chpair_dim1];
5014 /* ------------------------------ The end -------------------------------
5018 AdvApp2Var_SysBase::mgsomsg_("MMA2CFV", 7L);
5023 //=======================================================================
5024 //function : mma2ds1_
5026 //=======================================================================
5027 int AdvApp2Var_ApproxF2var::mma2ds1_(integer *ndimen,
5030 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
5045 /* System generated locals */
5046 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5047 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5048 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5049 fpntab_offset, i__1;
5051 /* Local variables */
5053 integer ibid1, ibid2, iuouv, nd;
5056 /* **********************************************************************
5061 /* Discretisation of function F(u,v) on the roots of Legendre polynoms. */
5065 /* FONCTION&,DISCRETISATION,&POINT */
5067 /* INPUT ARGUMENTS : */
5068 /* ------------------ */
5069 /* NDIMEN: Dimension of the space. */
5070 /* UINTFN: Limits of the interval of definition by u of the function */
5071 /* to be processed: (UINTFN(1),UINTFN(2)). */
5072 /* VINTFN: Limits of the interval of definition by v of the function */
5073 /* to be processed: (VINTFN(1),VINTFN(2)). */
5074 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5075 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5076 /* FONCNP is discretized by u. */
5077 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5078 /* FONCNP is discretized by v. */
5079 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5080 /* of Legendre of degree NBPNTU defined on (-1,1). */
5081 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5082 /* of Legendre of degree NBPNTV defined on (-1,1). */
5083 /* ISOFAV: Shows the type of iso of F(u,v) to be extracted to improve */
5084 /* the rapidity of calculation (has no influence on the form */
5086 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5087 /* with fixed u (with NBPNTV values different from v). */
5088 /* = 2, shows that it is necessaty to calculate the points of F(u,v) */
5089 /* with fixed v (with NBPNTU values different from u). */
5090 /* SOSOTB: Preinitialized table (input/output argument). */
5091 /* DISOTB: Preinitialized table (input/output argument). */
5092 /* SODITB: Preinitialized table (input/output argument). */
5093 /* DIDITB: Preinitialized table (input/output argument). */
5095 /* OUTPUT ARGUMENTS : */
5096 /* ------------------- */
5097 /* SOSOTB: Table where the terms */
5098 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5099 /* are added with ui and vj positive roots of Legendre polynom */
5100 /* of degree NBPNTU and NBPNTV respectively. */
5101 /* DISOTB: Table where the terms */
5102 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5103 /* are added with ui and vj positive roots of Legendre polynom */
5104 /* of degree NBPNTU and NBPNTV respectively. */
5105 /* SODITB: Table where the terms */
5106 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5107 /* are added with ui and vj positive roots of Legendre polynom */
5108 /* of degree NBPNTU and NBPNTV respectively. */
5109 /* DIDITB: Table where the terms */
5110 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5111 /* are added with ui and vj positive roots of Legendre polynom */
5112 /* of degree NBPNTU and NBPNTV respectively. */
5113 /* FPNTAB: Auxiliary table. */
5114 /* TTABLE: Auxiliary table. */
5115 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5116 /* the returned error code is equal to error code of FONCNP + 100. */
5118 /* COMMONS USED : */
5119 /* ---------------- */
5121 /* REFERENCES CALLED : */
5122 /* --------------------- */
5124 /* DESCRIPTION/NOTES/LIMITATIONS : */
5125 /* ----------------------------------- */
5126 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5127 /* where MA2FXK should be in the following form : */
5128 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,ISOFAV,TCONST,NBPTAB */
5129 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5130 /* with the following input arguments : */
5131 /* - NDIMEN is integer defined as the sum of dimensions of */
5132 /* sub-spaces (i.e. total dimension of the problem). */
5133 /* - UINTFN(2) is a table of 2 reals containing the interval */
5134 /* by u where the function to be approximated is defined */
5135 /* (so it is equal to UIFONC). */
5136 /* - VINTFN(2) is a table of 2 reals containing the interval */
5137 /* by v where the function to be approximated is defined */
5138 /* (so it is equal to VIFONC). */
5139 /* - ISOFAV, is 1 if it is necessary to calculate points with constant u, */
5140 /* is 2 if it is necessary to calculate points with constant v. */
5141 /* Any other value is an error. */
5142 /* - TCONST, real, value of the fixed parameter. Takes values */
5143 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5144 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5145 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5146 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5147 /* 'free' parameter of discretization (v if IISOFAV=1, */
5148 /* u if IISOFAV=2). */
5149 /* - IDERIU, integer, takes values between 0 (position) */
5150 /* and IORDRE(1) (partial derivative of the function by u */
5151 /* of order IORDRE(1) if IORDRE(1) > 0). */
5152 /* - IDERIV, integer, takes values between 0 (position) */
5153 /* and IORDRE(2) (partial derivative of the function by v */
5154 /* of order IORDRE(2) if IORDRE(2) > 0). */
5155 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5156 /* points of the derivative : */
5163 /* and the output arguments aret : */
5164 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5165 /* NBPTAB points calculated in FONCNP. */
5166 /* - IERCOD is, at output the error code of FONCNP. This code */
5167 /* (integer) should be strictly positive if there is a problem. */
5169 /* The input arguments SHOULD NOT be modified under FONCNP.
5172 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5173 /* values of UROOTB and VROOTB are consequently modified. */
5175 /* -->The results of discretisation are ranked in 4 tables */
5176 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5177 /* during the calculation of coefficients of the polynom of approximation. */
5179 /* When NBPNTU is uneven : */
5180 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5181 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5182 /* When NBPNTV is uneven : */
5183 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5184 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5185 /* When NBPNTU and NBPNTV are uneven : */
5186 /* term SOSOTB(0,0) contains F(0,0). */
5189 /* **********************************************************************
5191 /* Name of the routine */
5194 /* --------------------------- Initialization --------------------------
5197 /* Parameter adjustments */
5198 fpntab_dim1 = *ndimen;
5199 fpntab_offset = fpntab_dim1 + 1;
5200 fpntab -= fpntab_offset;
5204 diditb_dim1 = *nbpntu / 2 + 1;
5205 diditb_dim2 = *nbpntv / 2 + 1;
5206 diditb_offset = diditb_dim1 * diditb_dim2;
5207 diditb -= diditb_offset;
5208 soditb_dim1 = *nbpntu / 2;
5209 soditb_dim2 = *nbpntv / 2;
5210 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5211 soditb -= soditb_offset;
5212 disotb_dim1 = *nbpntu / 2;
5213 disotb_dim2 = *nbpntv / 2;
5214 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5215 disotb -= disotb_offset;
5216 sosotb_dim1 = *nbpntu / 2 + 1;
5217 sosotb_dim2 = *nbpntv / 2 + 1;
5218 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5219 sosotb -= sosotb_offset;
5224 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5226 AdvApp2Var_SysBase::mgenmsg_("MMA2DS1", 7L);
5229 if (*isofav < 1 || *isofav > 2) {
5235 /* **********************************************************************
5237 /* --------- Discretization by U on the roots of the polynom of ------ */
5238 /* --------------- Legendre of degree NBPNTU, iso-V by iso-V --------- */
5239 /* **********************************************************************
5243 mma2ds2_(ndimen, &uintfn[1], &vintfn[1], foncnp, nbpntu, nbpntv, &
5244 urootb[1], &vrootb[1], &iuouv, &sosotb[sosotb_offset], &
5245 disotb[disotb_offset], &soditb[soditb_offset], &diditb[
5246 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5248 /* ******************************************************************
5250 /* --------- Discretization by V on the roots of the polynom of ------ */
5251 /* --------------- Legendre of degree NBPNTV, iso-V by iso-V --------- */
5252 /* ******************************************************************
5256 /* --> Inversion of indices of tables */
5258 for (nd = 1; nd <= i__1; ++nd) {
5259 isz1 = *nbpntu / 2 + 1;
5260 isz2 = *nbpntv / 2 + 1;
5261 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5262 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5263 ibid1, &ibid2, iercod);
5267 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5268 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5269 ibid1, &ibid2, iercod);
5275 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5276 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5277 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5281 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5282 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5283 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5290 mma2ds2_(ndimen, &vintfn[1], &uintfn[1], foncnp, nbpntv, nbpntu, &
5291 vrootb[1], &urootb[1], &iuouv, &sosotb[sosotb_offset], &
5292 soditb[soditb_offset], &disotb[disotb_offset], &diditb[
5293 diditb_offset], &fpntab[fpntab_offset], &ttable[1], iercod);
5294 /* --> Inversion of indices of tables */
5296 for (nd = 1; nd <= i__1; ++nd) {
5297 isz1 = *nbpntv / 2 + 1;
5298 isz2 = *nbpntu / 2 + 1;
5299 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &isz1, &
5300 isz2, &isz2, &sosotb[nd * sosotb_dim2 * sosotb_dim1], &
5301 ibid1, &ibid2, iercod);
5305 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &diditb[nd * diditb_dim2 * diditb_dim1], &isz1, &
5306 isz2, &isz2, &diditb[nd * diditb_dim2 * diditb_dim1], &
5307 ibid1, &ibid2, iercod);
5313 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &soditb[(nd * soditb_dim2 + 1) * soditb_dim1 + 1],
5314 &isz1, &isz2, &isz2, &soditb[(nd * soditb_dim2 + 1) *
5315 soditb_dim1 + 1], &ibid1, &ibid2, iercod);
5319 AdvApp2Var_MathBase::mmfmtb1_(&isz1, &disotb[(nd * disotb_dim2 + 1) * disotb_dim1 + 1],
5320 &isz1, &isz2, &isz2, &disotb[(nd * disotb_dim2 + 1) *
5321 disotb_dim1 + 1], &ibid1, &ibid2, iercod);
5329 /* ------------------------------ The end -------------------------------
5335 AdvApp2Var_SysBase::maermsg_("MMA2DS1", iercod, 7L);
5338 AdvApp2Var_SysBase::mgsomsg_("MMA2DS1", 7L);
5343 //=======================================================================
5344 //function : mma2ds2_
5346 //=======================================================================
5347 int mma2ds2_(integer *ndimen,
5350 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
5366 /* System generated locals */
5367 integer sosotb_dim1, sosotb_dim2, sosotb_offset, disotb_dim1, disotb_dim2,
5368 disotb_offset, soditb_dim1, soditb_dim2, soditb_offset,
5369 diditb_dim1, diditb_dim2, diditb_offset, fpntab_dim1,
5370 fpntab_offset, i__1, i__2, i__3;
5372 /* Local variables */
5375 doublereal alinu, blinu, alinv, blinv, tcons;
5376 doublereal dbfn1[2], dbfn2[2];
5377 integer nuroo, nvroo, id, iu, iv;
5381 /* **********************************************************************
5386 /* Discretization of function F(u,v) on the roots of polynoms of Legendre. */
5390 /* FONCTION&,DISCRETISATION,&POINT */
5392 /* INPUT ARGUMENTS : */
5393 /* ------------------ */
5394 /* NDIMEN: Dimension of the space. */
5395 /* UINTFN: Limits of the interval of definition by u of the function */
5396 /* to be processed: (UINTFN(1),UINTFN(2)). */
5397 /* VINTFN: Limits of the interval of definition by v of the function */
5398 /* to be processed: (VINTFN(1),VINTFN(2)). */
5399 /* FONCNP: The NAME of the non-polynomial function to be processed. */
5400 /* NBPNTU: The degree of Legendre polynom on the roots which of */
5401 /* FONCNP is discretized by u. */
5402 /* NBPNTV: The degree of Legendre polynom on the roots which of */
5403 /* FONCNP is discretized by v. */
5404 /* UROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5405 /* of Legendre of degree NBPNTU defined on (-1,1). */
5406 /* VROOTB: Table of STRICTLY POSITIVE roots of the polynom */
5407 /* of Legendre of degree NBPNTV defined on (-1,1). */
5408 /* IIUOUV: Shows the type of iso of F(u,v) tom be extracted to improve the */
5409 /* rapidity of calculation (has no influence on the form of result) */
5410 /* = 1, shows that it is necessary to calculate the points of F(u,v) */
5411 /* with fixed u (so with NBPNTV values different from v). */
5412 /* = 2, shows that it is necessary to calculate the points of F(u,v) */
5413 /* with fixed v (so with NBPNTV values different from u). */
5414 /* SOSOTB: Preinitialized table (input/output argument). */
5415 /* DISOTB: Preinitialized table (input/output argument). */
5416 /* SODITB: Preinitialized table (input/output argument). */
5417 /* DIDITB: Preinitialized table (input/output argument). */
5419 /* OUTPUT ARGUMENTS : */
5420 /* ------------------- */
5421 /* SOSOTB: Table where the terms */
5422 /* F(ui,vj) + F(ui,-vj) + F(-ui,vj) + F(-ui,-vj) */
5423 /* are added with ui and vj positive roots of Legendre polynom */
5424 /* of degree NBPNTU and NBPNTV respectively. */
5425 /* DISOTB: Table where the terms */
5426 /* F(ui,vj) + F(ui,-vj) - F(-ui,vj) - F(-ui,-vj) */
5427 /* are added with ui and vj positive roots of Legendre polynom */
5428 /* of degree NBPNTU and NBPNTV respectively. */
5429 /* SODITB: Table where the terms */
5430 /* F(ui,vj) - F(ui,-vj) + F(-ui,vj) - F(-ui,-vj) */
5431 /* are added with ui and vj positive roots of Legendre polynom */
5432 /* of degree NBPNTU and NBPNTV respectively. */
5433 /* DIDITB: Table where the terms */
5434 /* F(ui,vj) - F(ui,-vj) - F(-ui,vj) + F(-ui,-vj) */
5435 /* are added with ui and vj positive roots of Legendre polynom */
5436 /* of degree NBPNTU and NBPNTV respectively. */
5437 /* FPNTAB: Auxiliary table. */
5438 /* TTABLE: Auxiliary table. */
5439 /* IERCOD: Error code >100 Pb in the evaluation of FONCNP, */
5440 /* the returned error code is equal to error code of FONCNP + 100. */
5442 /* COMMONS USED : */
5443 /* ---------------- */
5445 /* REFERENCES CALLED : */
5446 /* --------------------- */
5448 /* DESCRIPTION/NOTES/LIMITATIONS : */
5449 /* ----------------------------------- */
5450 /* --> The external function created by the caller of MA2F1K, MA2FDK */
5451 /* where MA2FXK should be in the following form : */
5452 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIIUOUV,TCONST,NBPTAB */
5453 /* ,TTABLE,IDERIU,IDERIV,PPNTAB,IERCOD) */
5454 /* with the following input arguments : */
5455 /* - NDIMEN is integer defined as the sum of dimensions of */
5456 /* sub-spaces (i.e. total dimension of the problem). */
5457 /* - UINTFN(2) is a table of 2 reals containing the interval */
5458 /* by u where the function to be approximated is defined */
5459 /* (so it is equal to UIFONC). */
5460 /* - VINTFN(2) is a table of 2 reals containing the interval */
5461 /* by v where the function to be approximated is defined */
5462 /* (so it is equal to VIFONC). */
5463 /* - IIIUOUV, is 1 if it is necessary to calculate points with constant u, */
5464 /* is 2 if it is necessary to calculate points with constant v. */
5465 /* Any other value is an error. */
5466 /* - TCONST, real, value of the fixed parameter. Takes values */
5467 /* in (UIFONC(1),UIFONC(2)) if ISOFAV = 1 or */
5468 /* ins (VIFONC(1),VIFONC(2)) if ISOFAV = 2. */
5469 /* - NBPTAB, integer. Shows the number of points to be calculated. */
5470 /* - TTABLE, a table of reals NBPTAB. These are the values of */
5471 /* 'free' parameter of discretization (v if IIIUOUV=1, */
5472 /* u if IIIUOUV=2). */
5473 /* - IDERIU, integer, takes values between 0 (position) */
5474 /* and IORDRE(1) (partial derivative of the function by u */
5475 /* of order IORDRE(1) if IORDRE(1) > 0). */
5476 /* - IDERIV, integer, takes values between 0 (position) */
5477 /* and IORDRE(2) (partial derivative of the function by v */
5478 /* of order IORDRE(2) if IORDRE(2) > 0). */
5479 /* If IDERIU=i and IDERIV=j, FONCNP should calculate the */
5480 /* points of the derivative : */
5487 /* and the output arguments aret : */
5488 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
5489 /* NBPTAB points calculated in FONCNP. */
5490 /* - IERCOD is, at output the error code of FONCNP. This code */
5491 /* (integer) should be strictly positive if there is a problem. */
5493 /* The input arguments SHOULD NOT be modified under FONCNP.
5496 /* -->As FONCNP is not forcedly defined in (-1,1)*(-1,1), the */
5497 /* values of UROOTB and VROOTB are consequently modified. */
5499 /* -->The results of discretisation are ranked in 4 tables */
5500 /* SOSOTB, DISOTB, SODITB and DIDITB to earn time */
5501 /* during the calculation of coefficients of the polynom of approximation. */
5503 /* When NBPNTU is uneven : */
5504 /* table SOSOTB(0,j) contains F(0,vj) + F(0,-vj), */
5505 /* table DIDITB(0,j) contains F(0,vj) - F(0,-vj), */
5506 /* When NBPNTV is uneven : */
5507 /* table SOSOTB(i,0) contains F(ui,0) + F(-ui,0), */
5508 /* table DIDITB(i,0) contains F(ui,0) - F(-ui,0), */
5509 /* When NBPNTU and NBPNTV are uneven : */
5510 /* term SOSOTB(0,0) contains F(0,0). */
5512 /* ATTENTION: These 4 tables are filled by varying the */
5513 /* 1st index first. So, the discretizations */
5514 /* of F(...,t) (for IIUOUV = 2) or of F(t,...) (IIUOUV = 1) */
5515 /* are stored in SOSOTB(...,t), SODITB(...,t), etc... */
5516 /* (this allows to gain important time). */
5517 /* It is required that the caller, in case of IIUOUV=1, */
5518 /* invert the roles of u and v, of SODITB and DISOTB BEFORE the */
5521 /* **********************************************************************
5524 /* Name of the routine */
5526 /* --> Indices of loops. */
5528 /* --------------------------- Initialization --------------------------
5531 /* Parameter adjustments */
5535 fpntab_dim1 = *ndimen;
5536 fpntab_offset = fpntab_dim1 + 1;
5537 fpntab -= fpntab_offset;
5539 diditb_dim1 = *nbpntu / 2 + 1;
5540 diditb_dim2 = *nbpntv / 2 + 1;
5541 diditb_offset = diditb_dim1 * diditb_dim2;
5542 diditb -= diditb_offset;
5543 soditb_dim1 = *nbpntu / 2;
5544 soditb_dim2 = *nbpntv / 2;
5545 soditb_offset = soditb_dim1 * (soditb_dim2 + 1) + 1;
5546 soditb -= soditb_offset;
5547 disotb_dim1 = *nbpntu / 2;
5548 disotb_dim2 = *nbpntv / 2;
5549 disotb_offset = disotb_dim1 * (disotb_dim2 + 1) + 1;
5550 disotb -= disotb_offset;
5551 sosotb_dim1 = *nbpntu / 2 + 1;
5552 sosotb_dim2 = *nbpntv / 2 + 1;
5553 sosotb_offset = sosotb_dim1 * sosotb_dim2;
5554 sosotb -= sosotb_offset;
5558 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5560 AdvApp2Var_SysBase::mgenmsg_("MMA2DS2", 7L);
5564 alinu = (uintfn[2] - uintfn[1]) / 2.;
5565 blinu = (uintfn[2] + uintfn[1]) / 2.;
5566 alinv = (vintfn[2] - vintfn[1]) / 2.;
5567 blinv = (vintfn[2] + vintfn[1]) / 2.;
5570 dbfn1[0] = vintfn[1];
5571 dbfn1[1] = vintfn[2];
5572 dbfn2[0] = uintfn[1];
5573 dbfn2[1] = uintfn[2];
5575 dbfn1[0] = uintfn[1];
5576 dbfn1[1] = uintfn[2];
5577 dbfn2[0] = vintfn[1];
5578 dbfn2[1] = vintfn[2];
5581 /* **********************************************************************
5583 /* -------- Discretization by U on the roots of Legendre polynom -------- */
5584 /* ---------------- of degree NBPNTU, with Vj fixed -------------------- */
5585 /* **********************************************************************
5588 nuroo = *nbpntu / 2;
5589 nvroo = *nbpntv / 2;
5590 jdec = (*nbpntu + 1) / 2;
5592 /* ----------- Loading of parameters of discretization by U ------------- */
5595 for (iu = 1; iu <= i__1; ++iu) {
5596 ttable[iu] = blinu + alinu * urootb[iu];
5600 /* -------------- For Vj fixed, negative root of Legendre ------------- */
5603 for (iv = 1; iv <= i__1; ++iv) {
5604 tcons = blinv + alinv * vrootb[iv];
5605 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5606 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5607 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5612 for (id = 1; id <= i__2; ++id) {
5614 for (iu = 1; iu <= i__3; ++iu) {
5615 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5616 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5617 sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1]
5618 = sosotb[iu + (nvroo - iv + 1 + id * sosotb_dim2) *
5619 sosotb_dim1] + up + um;
5620 disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) * disotb_dim1]
5621 = disotb[iu + (nvroo - iv + 1 + id * disotb_dim2) *
5622 disotb_dim1] + up - um;
5623 soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) * soditb_dim1]
5624 = soditb[iu + (nvroo - iv + 1 + id * soditb_dim2) *
5625 soditb_dim1] - up - um;
5626 diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1]
5627 = diditb[iu + (nvroo - iv + 1 + id * diditb_dim2) *
5628 diditb_dim1] - up + um;
5631 if (*nbpntu % 2 != 0) {
5632 up = fpntab[id + jdec * fpntab_dim1];
5633 sosotb[(nvroo - iv + 1 + id * sosotb_dim2) * sosotb_dim1] +=
5635 diditb[(nvroo - iv + 1 + id * diditb_dim2) * diditb_dim1] -=
5643 /* --------- For Vj = 0 (uneven NBPNTV), discretization by U ----------- */
5645 if (*nbpntv % 2 != 0) {
5647 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5648 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5649 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5654 for (id = 1; id <= i__1; ++id) {
5656 for (iu = 1; iu <= i__2; ++iu) {
5657 up = fpntab[id + (jdec + iu) * fpntab_dim1];
5658 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5659 sosotb[iu + id * sosotb_dim2 * sosotb_dim1] = sosotb[iu + id *
5660 sosotb_dim2 * sosotb_dim1] + up + um;
5661 diditb[iu + id * diditb_dim2 * diditb_dim1] = diditb[iu + id *
5662 diditb_dim2 * diditb_dim1] + up - um;
5665 if (*nbpntu % 2 != 0) {
5666 up = fpntab[id + jdec * fpntab_dim1];
5667 sosotb[id * sosotb_dim2 * sosotb_dim1] += up;
5673 /* -------------- For Vj fixed, positive root of Legendre ------------- */
5676 for (iv = 1; iv <= i__1; ++iv) {
5677 tcons = alinv * vrootb[(*nbpntv + 1) / 2 + iv] + blinv;
5678 (*const_cast <AdvApp2Var_EvaluatorFunc2Var*> (&foncnp)).Evaluate (
5679 ndimen, dbfn1, dbfn2, iiuouv, &tcons, nbpntu, &
5680 ttable[1], &c__0, &c__0, &fpntab[fpntab_offset], iercod);
5685 for (id = 1; id <= i__2; ++id) {
5687 for (iu = 1; iu <= i__3; ++iu) {
5688 up = fpntab[id + (iu + jdec) * fpntab_dim1];
5689 um = fpntab[id + (nuroo - iu + 1) * fpntab_dim1];
5690 sosotb[iu + (iv + id * sosotb_dim2) * sosotb_dim1] = sosotb[
5691 iu + (iv + id * sosotb_dim2) * sosotb_dim1] + up + um;
5692 disotb[iu + (iv + id * disotb_dim2) * disotb_dim1] = disotb[
5693 iu + (iv + id * disotb_dim2) * disotb_dim1] + up - um;
5694 soditb[iu + (iv + id * soditb_dim2) * soditb_dim1] = soditb[
5695 iu + (iv + id * soditb_dim2) * soditb_dim1] + up + um;
5696 diditb[iu + (iv + id * diditb_dim2) * diditb_dim1] = diditb[
5697 iu + (iv + id * diditb_dim2) * diditb_dim1] + up - um;
5700 if (*nbpntu % 2 != 0) {
5701 up = fpntab[id + jdec * fpntab_dim1];
5702 sosotb[(iv + id * sosotb_dim2) * sosotb_dim1] += up;
5703 diditb[(iv + id * diditb_dim2) * diditb_dim1] += up;
5710 /* ------------------------------ The end -------------------------------
5716 AdvApp2Var_SysBase::maermsg_("MMA2DS2", iercod, 7L);
5719 AdvApp2Var_SysBase::mgsomsg_("MMA2DS2", 7L);
5724 //=======================================================================
5725 //function : mma2er1_
5727 //=======================================================================
5728 int mma2er1_(integer *ndjacu,
5744 /* System generated locals */
5745 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
5748 /* Local variables */
5753 doublereal bid0, bid1;
5755 /* **********************************************************************
5760 /* Calculate max approximation error done when */
5761 /* the coefficients of PATJAC such that the degree by U varies between */
5762 /* MINDGU and MAXDGU and the degree by V varies between MINDGV and MAXDGV are removed. */
5766 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5768 /* INPUT ARGUMENTS : */
5769 /* ------------------ */
5770 /* NDJACU: Dimension by U of table PATJAC. */
5771 /* NDJACV: Dimension by V of table PATJAC. */
5772 /* NDIMEN: Dimension of the space. */
5773 /* MINDGU: Lower limit of index by U of coeff. of PATJAC to be taken into account. */
5774 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5775 /* MINDGV: Lower limit of index by V of coeff. of PATJAC to be taken into account. */
5776 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5777 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5778 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5779 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5780 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5781 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5782 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5783 /* PATJAC: Table of coeff. of square of approximation with */
5784 /* constraints of order IORDRU by U and IORDRV by V. */
5785 /* VECERR: Auxiliary vector. */
5786 /* ERREUR: MAX Error commited during removal of ALREADY CALCULATED coeff of PATJAC */
5788 /* OUTPUT ARGUMENTS : */
5789 /* ------------------- */
5790 /* ERREUR: MAX Error commited during removal of coeff of PATJAC */
5791 /* of indices from MINDGU to MAXDGU by U and from MINDGV to MAXDGV by V */
5792 /* THEN the already calculated error. */
5794 /* COMMONS USED : */
5795 /* ---------------- */
5797 /* REFERENCES CALLED : */
5798 /* --------------------- */
5800 /* DESCRIPTION/NOTES/LIMITATIONS : */
5801 /* ----------------------------------- */
5802 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5803 /* approximation of F(U,V). The indices i and j show the degree */
5804 /* by U and by V of base polynoms. These polynoms have the form: */
5806 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5808 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5809 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5811 /* The contribution to the error of term Cij when it is */
5812 /* removed from PATJAC is increased by: */
5814 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5816 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5818 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5822 /* ***********************************************************************
5824 /* Name of the routine */
5827 /* ----------------------------- Initialisations ------------------------
5830 /* Parameter adjustments */
5832 patjac_dim1 = *ndjacu + 1;
5833 patjac_dim2 = *ndjacv + 1;
5834 patjac_offset = patjac_dim1 * patjac_dim2;
5835 patjac -= patjac_offset;
5838 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
5840 AdvApp2Var_SysBase::mgenmsg_("MMA2ER1", 7L);
5843 minu = (*iordru + 1) << 1;
5844 minv = (*iordrv + 1) << 1;
5846 /* ------------------- Calculate the increment of the max error --------------- */
5847 /* ----- during the removal of the coeffs of indices from MINDGU to MAXDGU ---- */
5848 /* ---------------- by U and indices from MINDGV to MAXDGV by V --------------- */
5851 for (nd = 1; nd <= i__1; ++nd) {
5854 for (jj = *mindgv; jj <= i__2; ++jj) {
5857 for (ii = *mindgu; ii <= i__3; ++ii) {
5858 bid0 += (d__1 = patjac[ii + (jj + nd * patjac_dim2) *
5859 patjac_dim1], advapp_abs(d__1)) * xmaxju[ii - minu];
5862 bid1 = bid0 * xmaxjv[jj - minv] + bid1;
5870 /* ----------------------- Calculate the max error ----------------------*/
5872 bid1 = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
5876 *erreur = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
5878 /* ------------------------- The end ------------------------------------
5882 AdvApp2Var_SysBase::mgsomsg_("MMA2ER1", 7L);
5887 //=======================================================================
5888 //function : mma2er2_
5890 //=======================================================================
5891 int mma2er2_(integer *ndjacu,
5903 doublereal *epmscut,
5910 /* System generated locals */
5911 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2;
5914 /* Local variables */
5917 integer i2rdu, i2rdv;
5918 doublereal errnu, errnv;
5919 integer ii, nd, jj, nu, nv;
5920 doublereal bid0, bid1;
5922 /* **********************************************************************
5927 /* Remove coefficients of PATJAC to obtain the minimum degree */
5928 /* by U and V checking the imposed tolerance. */
5932 /* TOUS,AB_SPECIFI:: CARREAU&,CALCUL,&ERREUR */
5934 /* INPUT ARGUMENTS : */
5935 /* ------------------ */
5936 /* NDJACU: Degree by U of table PATJAC. */
5937 /* NDJACV: Degree by V of table PATJAC. */
5938 /* NDIMEN: Dimension of the space. */
5939 /* MINDGU: Limit of index by U of coeff. of PATJAC to be PRESERVED (should be >=0). */
5940 /* MAXDGU: Upper limit of index by U of coeff. of PATJAC to be taken into account. */
5941 /* MINDGV: Limit of index by V of coeff. of PATJAC to be PRESERVED (should be >=0). */
5942 /* MAXDGV: Upper limit of index by V of coeff. of PATJAC to be taken into account. */
5943 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5944 /* IORDRV: Order of continuity by U provided by square PATJAC (from -1 to 2) */
5945 /* XMAXJU: Maximum value of Jacobi polynoms of order IORDRU, */
5946 /* from degree 0 to MAXDGU - 2*(IORDU+1) */
5947 /* XMAXJV: Maximum value of Jacobi polynoms of order IORDRV, */
5948 /* from degree 0 to MAXDGV - 2*(IORDV+1) */
5949 /* PATJAC: Table of coeff. of square of approximation with */
5950 /* constraints of order IORDRU by U and IORDRV by V. */
5951 /* EPMSCUT: Tolerance of approximation. */
5952 /* VECERR: Auxiliary vector. */
5953 /* ERREUR: MAX Error commited ALREADY CALCULATED */
5955 /* OUTPUT ARGUMENTS : */
5956 /* ------------------- */
5957 /* ERREUR: MAX Error commited by preserving only coeff of PATJAC */
5958 /* of indices from 0 to NEWDGU by U and from 0 to NEWDGV by V */
5959 /* PLUS the already calculated error. */
5960 /* NEWDGU: Min. Degree by U such as the square of approximation */
5961 /* could check the tolerance. There is always NEWDGU >= MINDGU >= 0. */
5962 /* NEWDGV: Min. Degree by V such as the square of approximation */
5963 /* could check the tolerance. There is always NEWDGV >= MINDGV >= 0. */
5966 /* COMMONS USED : */
5967 /* ---------------- */
5969 /* REFERENCES CALLED : */
5970 /* --------------------- */
5972 /* DESCRIPTION/NOTES/LIMITATIONS : */
5973 /* ----------------------------------- */
5974 /* Table PATJAC is the place of storage of coeff. Cij of the square of */
5975 /* approximation of F(U,V). The indices i and j show the degree */
5976 /* by U and by V of base polynoms. These polynoms have the form: */
5978 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
5980 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
5981 /* IORDRU+1 (the same by V by replacing U u V in the expression above). */
5983 /* The contribution to the error of term Cij when it is */
5984 /* removed from PATJAC is increased by: */
5986 /* DABS(Cij)*XMAXJU(i-2*(IORDRU+1))*XMAXJV(J-2*(IORDRV+1)) where */
5988 /* XMAXJU(i-2*(IORDRU+1) = ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U),
5990 /* XMAXJV(i-2*(IORDRV+1) = ((1 - V*V)**(IORDRV+1)).J(j-2*(IORDRV+1)(V).
5994 /* **********************************************************************
5996 /* Name of the routine */
5999 /* ----------------------------- Initialisations ------------------------
6002 /* Parameter adjustments */
6004 patjac_dim1 = *ndjacu + 1;
6005 patjac_dim2 = *ndjacv + 1;
6006 patjac_offset = patjac_dim1 * patjac_dim2;
6007 patjac -= patjac_offset;
6010 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
6012 AdvApp2Var_SysBase::mgenmsg_("MMA2ER2", 7L);
6015 i2rdu = (*iordru + 1) << 1;
6016 i2rdv = (*iordrv + 1) << 1;
6020 /* **********************************************************************
6022 /* -------------------- Cutting of oefficients ------------------------
6024 /* **********************************************************************
6029 /* ------------------- Calculate the increment of max error --------------- */
6030 /* ----- during the removal of coeff. of indices from MINDGU to MAXDGU ------ */
6031 /* ---------------- by U, the degree by V is fixed to NV -----------------
6036 bid0 = xmaxjv[nv - i2rdv];
6038 for (nd = 1; nd <= i__1; ++nd) {
6041 for (ii = i2rdu; ii <= i__2; ++ii) {
6042 bid1 += (d__1 = patjac[ii + (nv + nd * patjac_dim2) *
6043 patjac_dim1], advapp_abs(d__1)) * xmaxju[ii - i2rdu] * bid0;
6050 vecerr[1] = *epmscut * 2;
6052 errnv = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6054 /* ------------------- Calculate the increment of max error --------------- */
6055 /* ----- during the removal of coeff. of indices from MINDGV to MAXDGV ------ */
6056 /* ---------------- by V, the degree by U is fixed to NU -----------------
6061 bid0 = xmaxju[nu - i2rdu];
6063 for (nd = 1; nd <= i__1; ++nd) {
6066 for (jj = i2rdv; jj <= i__2; ++jj) {
6067 bid1 += (d__1 = patjac[nu + (jj + nd * patjac_dim2) *
6068 patjac_dim1], advapp_abs(d__1)) * xmaxjv[jj - i2rdv] * bid0;
6075 vecerr[1] = *epmscut * 2;
6077 errnu = AdvApp2Var_MathBase::mzsnorm_(ndimen, &vecerr[1]);
6079 /* ----------------------- Calculate the max error ----------------------
6085 errnu = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6087 errnv = AdvApp2Var_MathBase::mzsnorm_(&nd, vaux);
6089 if (errnu > errnv) {
6090 if (errnv < *epmscut) {
6097 if (errnu < *epmscut) {
6107 /* -------------------------- Return the degrees -------------------
6111 *newdgu = advapp_max(nu,1);
6112 *newdgv = advapp_max(nv,1);
6114 /* ----------------------------------- The end --------------------------
6118 AdvApp2Var_SysBase::mgsomsg_("MMA2ER2", 7L);
6123 //=======================================================================
6124 //function : mma2fnc_
6126 //=======================================================================
6127 int AdvApp2Var_ApproxF2var::mma2fnc_(integer *ndimen,
6131 const AdvApp2Var_EvaluatorFunc2Var& foncnp,
6157 /* System generated locals */
6158 integer courbe_dim1, courbe_dim2, courbe_offset, somtab_dim1, somtab_dim2,
6159 somtab_offset, diftab_dim1, diftab_dim2, diftab_offset,
6160 contr1_dim1, contr1_dim2, contr1_offset, contr2_dim1, contr2_dim2,
6161 contr2_offset, errmax_dim1, errmax_offset, errmoy_dim1,
6162 errmoy_offset, i__1;
6165 /* Local variables */
6168 integer ideb1, ibid1, ibid2, ncfja, ndgre, ilong,
6170 doublereal* wrkar = 0;
6173 doublereal uvpav[4] /* was [2][2] */;
6177 doublereal uv11[4] /* was [2][2] */;
6180 integer isz1, isz2, isz3, isz4, isz5;
6181 intptr_t ipt1, ipt2, ipt3, ipt4, ipt5,iptt, jptt;
6183 /* **********************************************************************
6188 /* Approximation of a limit of non polynomial function F(u,v) */
6189 /* (in the space of dimension NDIMEN) by SEVERAL */
6190 /* polynomial curves, by the method of least squares. The parameter of the function is preserved. */
6194 /* TOUS, AB_SPECIFI :: FONCTION&,EXTREMITE&, APPROXIMATION, &COURBE. */
6196 /* INPUT ARGUMENTS : */
6197 /* ----------------- */
6198 /* NDIMEN: Total Dimension of the space (sum of dimensions */
6199 /* of sub-spaces) */
6200 /* NBSESP: Number of "independent" sub-spaces. */
6201 /* NDIMSE: Table of dimensions of sub-spaces. */
6202 /* UVFONC: Limits of the interval (a,b)x(c,d) of definition of the */
6203 /* function to be approached by U (UVFONC(*,1) contains (a,b)) */
6204 /* and by V (UVFONC(*,2) contains (c,d)). */
6205 /* FONCNP: External function of position on the non polynomial function to be approached. */
6206 /* TCONST: Value of isoparameter of F(u,v) to be discretized. */
6207 /* ISOFAV: Type of chosen iso, = 1, shose that discretization is with u */
6208 /* fixed; = 2, shows that v is fixed. */
6209 /* NBROOT: Nb of points of discretisation of the iso, extremities not included. */
6210 /* ROOTLG: Table of roots of the polynom of Legendre defined on */
6211 /* (-1,1), of degree NBROOT. */
6212 /* IORDRE: Order of constraint at the extremities of the limit */
6213 /* -1 = no constraints, */
6214 /* 0 = constraints of passage to limits (i.e. C0), */
6215 /* 1 = C0 + constraints of 1st derivatives (i.e. C1), */
6216 /* 2 = C1 + constraints of 2nd derivatives (i.e. C2). */
6217 /* IDERIV: Order of derivative of the limit. */
6218 /* NDGJAC: Degree of serial development to be used for calculation in */
6219 /* the Jacobi base. */
6220 /* NBCRMX: Max Nb of curves to be created. */
6221 /* NCFLIM: Max Nb of coeff of the polynomial curve */
6222 /* of approximation (should be above or equal to */
6223 /* 2*IORDRE+2 and below or equal to 50). */
6224 /* EPSAPR: Table of required errors of approximation */
6225 /* sub-space by sub-space. */
6227 /* OUTPUT ARGUMENTS : */
6228 /* ------------------- */
6229 /* NCOEFF: Number of significative coeff of calculated curves. */
6230 /* COURBE: Table of coeff. of calculated polynomial curves. */
6231 /* Should be dimensioned in (NCFLIM,NDIMEN,NBCRMX). */
6232 /* These curves are ALWAYS parametrized in (-1,1). */
6233 /* NBCRBE: Nb of calculated curves. */
6234 /* SOMTAB: For F defined on (-1,1) (otherwise rescale the */
6235 /* parameters), this is the table of sums F(u,vj) + F(u,-vj)
6237 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6238 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6239 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6240 /* v of order IDERIV are taken). */
6241 /* DIFTAB: For F defined on (-1,1) (otherwise rescale the */
6242 /* parameters), this is the table of sums F(u,vj) - F(u,-vj)
6244 /* if ISOFAV = 1 (and IDERIV=0, otherwise the derivatives */
6245 /* by u of order IDERIV are taken) or sumes F(ui,v) + F(-ui,v) if */
6246 /* ISOFAV = 2 (and IDERIV=0, otherwise the derivatives by */
6247 /* v of order IDERIV are taken). */
6248 /* CONTR1: Contains the coordinates of the left extremity of the iso */
6249 /* and of its derivatives till order IORDRE */
6250 /* CONTR2: Contains the coordinates of the right extremity of the iso */
6251 /* and of its derivatives till order IORDRE */
6252 /* TABDEC: Table of NBCRBE+1 parameters of cut of UVFONC(1:2,1)
6254 /* if ISOFAV=2, or of UVFONC(1:2,2) if ISOFAV=1. */
6255 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6256 /* committed in the approximation of FONCNP by NBCRBE curves. */
6257 /* ERRMOY: Table of AVERAGE errors (sub-space by sub-space) */
6258 /* committed in the approximation of FONCNP by NBCRBE curves. */
6259 /* IERCOD: Error code: */
6260 /* -1 = ERRMAX > EPSAPR for at least one sub-space. */
6261 /* (the resulting curves of at least mathematic degree NCFLIM-1 */
6262 /* are calculated). */
6263 /* 0 = Everything is ok. */
6264 /* 1 = Pb of incoherence of inputs. */
6265 /* 10 = Pb of calculation of the interpolation of constraints. */
6266 /* 13 = Pb in the dynamic allocation. */
6267 /* 33 = Pb in the data recuperation from block data */
6268 /* of coeff. of integration by GAUSS method. */
6269 /* >100 Pb in the evaluation of FONCNP, the returned error code */
6270 /* is equal to the error code of FONCNP + 100. */
6272 /* COMMONS USED : */
6273 /* ---------------- */
6275 /* REFERENCES CALLED : */
6276 /* ----------------------- */
6278 /* DESCRIPTION/NOTES/LIMITATIONS : */
6279 /* ----------------------------------- */
6280 /* --> The approximation part is done in the space of dimension */
6281 /* NDIMEN (the sum of dimensions of sub-spaces). For example : */
6282 /* If NBSESP=2 and NDIMSE(1)=3, NDIMSE(2)=2, there is smoothing with */
6283 /* NDIMEN=5. The result (in COURBE(NDIMEN,NCOEFF,i) ), will be */
6284 /* composed of the result of smoothing of 3D function in */
6285 /* COURBE(1:3,1:NCOEFF,i) and of smoothing of 2D function in */
6286 /* COURBE(4:5,1:NCOEFF,i). */
6288 /* --> Routine FONCNP should be declared EXTERNAL in the program */
6289 /* calling MMA2FNC. */
6291 /* --> Function FONCNP, declared externally, should be declared */
6292 /* IMPERATIVELY in form : */
6293 /* SUBROUTINE FONCNP(NDIMEN,UINTFN,VINTFN,IIUOUV,TCONST,NBPTAB */
6294 /* ,TTABLE,IDERIU,IDERIV,IERCOD) */
6295 /* where the input arguments are : */
6296 /* - NDIMEN is integer defined as the sum of dimensions of */
6297 /* sub-spaces (i.e. total dimension of the problem). */
6298 /* - UINTFN(2) is a table of 2 reals containing the interval */
6299 /* by u where the function to be approximated is defined */
6300 /* (so it is equal to UIFONC). */
6301 /* - VINTFN(2) is a table of 2 reals containing the interval */
6302 /* by v where the function to be approximated is defined */
6303 /* (so it is equal to VIFONC). */
6304 /* - IIUOUV, shows that the points to be calculated have a constant U */
6305 /* (IIUOUV=1) or a constant V (IIUOUV=2). */
6306 /* - TCONST, real, value of the fixed discretisation parameter. Takes values */
6307 /* in (UINTFN(1),UINTFN(2)) if IIUOUV=1, */
6308 /* or in (VINTFN(1),VINTFN(2)) if IIUOUV=2. */
6309 /* - NBPTAB, the nb of point of discretisation following the free variable */
6310 /* : V if IIUOUV=1 or U if IIUOUV = 2. */
6311 /* - TTABLE, Table of NBPTAB parametres of discretisation. . */
6312 /* - IDERIU, integer, takes values between 0 (position) */
6313 /* and IORDREU (partial derivative of the function by u */
6314 /* of order IORDREU if IORDREU > 0). */
6315 /* - IDERIV, integer, takes values between 0 (position) */
6316 /* and IORDREV (partial derivative of the function by v */
6317 /* of order IORDREV if IORDREV > 0). */
6318 /* and the output arguments are : */
6319 /* - FPNTAB(NDIMEN,NBPTAB) contains, at output, the table of */
6320 /* NBPTAB points calculated in FONCNP. */
6321 /* - IERCOD is, at output the error code of FONCNP. This code */
6322 /* (integer) should be strictly positive if there is a problem. */
6324 /* The input arguments SHOULD NOT BE modified under FONCNP.
6327 /* --> If IERCOD=-1, the required precision can't be reached (ERRMAX */
6328 /* is above EPSAPR on at least one sub-space), but
6330 /* one gives the best possible result for NCFLIM and EPSAPR */
6331 /* chosen by the user. In this case (and for IERCOD=0), there is a solution. */
6334 /* **********************************************************************
6336 /* Name of the routine */
6338 /* Parameter adjustments */
6343 errmoy_dim1 = *nbsesp;
6344 errmoy_offset = errmoy_dim1 + 1;
6345 errmoy -= errmoy_offset;
6346 errmax_dim1 = *nbsesp;
6347 errmax_offset = errmax_dim1 + 1;
6348 errmax -= errmax_offset;
6349 contr2_dim1 = *ndimen;
6350 contr2_dim2 = *iordre + 2;
6351 contr2_offset = contr2_dim1 * (contr2_dim2 + 1) + 1;
6352 contr2 -= contr2_offset;
6353 contr1_dim1 = *ndimen;
6354 contr1_dim2 = *iordre + 2;
6355 contr1_offset = contr1_dim1 * (contr1_dim2 + 1) + 1;
6356 contr1 -= contr1_offset;
6357 diftab_dim1 = *nbroot / 2 + 1;
6358 diftab_dim2 = *ndimen;
6359 diftab_offset = diftab_dim1 * (diftab_dim2 + 1);
6360 diftab -= diftab_offset;
6361 somtab_dim1 = *nbroot / 2 + 1;
6362 somtab_dim2 = *ndimen;
6363 somtab_offset = somtab_dim1 * (somtab_dim2 + 1);
6364 somtab -= somtab_offset;
6366 courbe_dim1 = *ncflim;
6367 courbe_dim2 = *ndimen;
6368 courbe_offset = courbe_dim1 * (courbe_dim2 + 1) + 1;
6369 courbe -= courbe_offset;
6370 AdvApp2Var_SysBase anAdvApp2Var_SysBase;
6373 ibb = AdvApp2Var_SysBase::mnfndeb_();
6375 AdvApp2Var_SysBase::mgenmsg_("MMA2FNC", 7L);
6380 /* ---------------- Set to zero the coefficients of CURVE --------------
6383 ilong = *ndimen * *ncflim * *nbcrmx;
6384 AdvApp2Var_SysBase::mvriraz_(&ilong, &courbe[courbe_offset]);
6386 /* **********************************************************************
6388 /* -------------------------- Checking of entries ------------------
6390 /* **********************************************************************
6393 AdvApp2Var_MathBase::mmveps3_(&eps3);
6394 if ((d__1 = uvfonc[4] - uvfonc[3], advapp_abs(d__1)) < eps3) {
6397 if ((d__1 = uvfonc[6] - uvfonc[5], advapp_abs(d__1)) < eps3) {
6406 /* ********************************************************************** */
6407 /* ------------- Preparation of parameters of discretisation ----------- */
6408 /* **********************************************************************
6411 /* -- Allocation of a table of parameters and points of discretisation -- */
6412 /* --> For the parameters of discretisation. */
6414 /* --> For the points of discretisation in MMA1FDI and MMA1CDI and
6416 /* the auxiliary curve for MMAPCMP */
6417 ibid1 = *ndimen * (*nbroot + 2);
6418 ibid2 = ((*iordre + 1) << 1) * *nbroot;
6419 isz2 = advapp_max(ibid1,ibid2);
6420 ibid1 = (((*ncflim - 1) / 2 + 1) << 1) * *ndimen;
6421 isz2 = advapp_max(ibid1,isz2);
6422 /* --> To return the polynoms of hermit. */
6423 isz3 = ((*iordre + 1) << 2) * (*iordre + 1);
6424 /* --> For the Gauss coeff. of integration. */
6425 isz4 = (*nbroot / 2 + 1) * (*ndgjac + 1 - ((*iordre + 1) << 1));
6426 /* --> For the coeff of the curve in the base of Jacobi */
6427 isz5 = (*ndgjac + 1) * *ndimen;
6429 ndwrk = isz1 + isz2 + isz3 + isz4 + isz5;
6430 anAdvApp2Var_SysBase.mcrrqst_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6433 /* --> For the parameters of discretisation (NBROOT+2 extremities). */
6435 /* --> For the points of discretisation FPNTAB(NDIMEN,NBROOT+2), */
6436 /* FPNTAB(NBROOT,2*(IORDRE+1)) and for WRKAR of MMAPCMP. */
6438 /* --> For the polynoms of Hermit */
6440 /* --> For the Gauss coeff of integration. */
6442 /* --> For the curve in Jacobi. */
6445 /* ------------------ Initialisation of management of cuts ---------
6449 uvpav[0] = uvfonc[3];
6450 uvpav[1] = uvfonc[4];
6451 tabdec[0] = uvfonc[5];
6452 tabdec[1] = uvfonc[6];
6453 } else if (*isofav == 2) {
6454 tabdec[0] = uvfonc[3];
6455 tabdec[1] = uvfonc[4];
6456 uvpav[2] = uvfonc[5];
6457 uvpav[3] = uvfonc[6];
6465 /* **********************************************************************
6467 /* APPROXIMATION WITH CUTS */
6468 /* **********************************************************************
6472 /* --> When the top is reached, this is the end ! */
6473 if (nupil - *nbcrbe == 0) {
6478 uvpav[2] = tabdec[*nbcrbe];
6479 uvpav[3] = tabdec[*nbcrbe + 1];
6480 } else if (*isofav == 2) {
6481 uvpav[0] = tabdec[*nbcrbe];
6482 uvpav[1] = tabdec[*nbcrbe + 1];
6487 /* -------------------- Normalization of parameters -------------------- */
6489 mma1nop_(nbroot, &rootlg[1], uvpav, isofav, &wrkar[ipt1], &ier);
6494 /* -------------------- Discretisation of FONCNP ------------------------ */
6496 mma1fdi_(ndimen, uvpav, foncnp, isofav, tconst, nbroot, &wrkar[ipt1],
6497 iordre, ideriv, &wrkar[ipt2], &somtab[(ncb1 * somtab_dim2 + 1) *
6498 somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6499 contr1[(ncb1 * contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1
6500 * contr2_dim2 + 1) * contr2_dim1 + 1], iercod);
6505 /* -----------Cut the discretisation of constraints ------------*/
6508 mma1cdi_(ndimen, nbroot, &rootlg[1], iordre, &contr1[(ncb1 *
6509 contr1_dim2 + 1) * contr1_dim1 + 1], &contr2[(ncb1 *
6510 contr2_dim2 + 1) * contr2_dim1 + 1], &somtab[(ncb1 *
6511 somtab_dim2 + 1) * somtab_dim1], &diftab[(ncb1 * diftab_dim2
6512 + 1) * diftab_dim1], &wrkar[ipt2], &wrkar[ipt3], &ier);
6518 /* **********************************************************************
6520 /* -------------------- Calculate the curve of approximation -------------
6522 /* **********************************************************************
6525 mma1jak_(ndimen, nbroot, iordre, ndgjac, &somtab[(ncb1 * somtab_dim2 + 1)
6526 * somtab_dim1], &diftab[(ncb1 * diftab_dim2 + 1) * diftab_dim1], &
6527 wrkar[ipt4], &wrkar[ipt5], &ier);
6532 /* **********************************************************************
6534 /* ---------------- Add polynom of interpolation -------------------
6536 /* **********************************************************************
6540 mma1cnt_(ndimen, iordre, &contr1[(ncb1 * contr1_dim2 + 1) *
6541 contr1_dim1 + 1], &contr2[(ncb1 * contr2_dim2 + 1) *
6542 contr2_dim1 + 1], &wrkar[ipt3], ndgjac, &wrkar[ipt5]);
6545 /* **********************************************************************
6547 /* --------------- Calculate Max and Average error ----------------------
6549 /* **********************************************************************
6552 mma1fer_(ndimen, nbsesp, &ndimse[1], iordre, ndgjac, &wrkar[ipt5], ncflim,
6553 &epsapr[1], &wrkar[ipt2], &errmax[ncb1 * errmax_dim1 + 1], &
6554 errmoy[ncb1 * errmoy_dim1 + 1], &ncoeff[ncb1], &ier);
6559 if (ier == 0 || (ier == -1 && nupil == *nbcrmx)) {
6561 /* ******************************************************************
6563 /* ----------------------- Compression du resultat ------------------
6565 /* ******************************************************************
6571 ncfja = *ndgjac + 1;
6572 /* -> Compression of result in WRKAR(IPT2) */
6575 AdvApp2Var_MathBase::mmapcmp_(ndimen,
6576 &ncfja, &ncoeff[ncb1], &wrkar[ipt5], &wrkar[ipt2]);
6578 AdvApp2Var_MathBase::mmapcmp_((integer*)ndimen,
6584 ilong = *ndimen * *ncflim;
6585 AdvApp2Var_SysBase::mvriraz_(&ilong, &wrkar[ipt5]);
6586 /* -> Passage to canonic base (-1,1) (result in WRKAR(IPT5)).
6588 ndgre = ncoeff[ncb1] - 1;
6590 for (nd = 1; nd <= i__1; ++nd) {
6591 iptt = ipt2 + ((nd - 1) << 1) * (ndgre / 2 + 1);
6592 jptt = ipt5 + (nd - 1) * ncoeff[ncb1];
6593 AdvApp2Var_MathBase::mmjacan_(iordre, &ndgre, &wrkar[iptt], &wrkar[jptt]);
6597 /* -> Store the calculated curve */
6599 AdvApp2Var_MathBase::mmfmca8_(&ncoeff[ncb1], ndimen, &ibid1, ncflim, ndimen, &ibid1, &
6600 wrkar[ipt5], &courbe[(ncb1 * courbe_dim2 + 1) * courbe_dim1 +
6603 /* -> Before normalization of constraints on (-1,1), recalculate */
6604 /* the true constraints. */
6606 for (ii = 0; ii <= i__1; ++ii) {
6607 mma1noc_(uv11, ndimen, &ii, &contr1[(ii + 1 + ncb1 * contr1_dim2)
6608 * contr1_dim1 + 1], uvpav, isofav, ideriv, &contr1[(ii +
6609 1 + ncb1 * contr1_dim2) * contr1_dim1 + 1]);
6610 mma1noc_(uv11, ndimen, &ii, &contr2[(ii + 1 + ncb1 * contr2_dim2)
6611 * contr2_dim1 + 1], uvpav, isofav, ideriv, &contr2[(ii +
6612 1 + ncb1 * contr2_dim2) * contr2_dim1 + 1]);
6616 ibid1 = (*nbroot / 2 + 1) * *ndimen;
6617 mma1noc_(uv11, &ibid1, &ii, &somtab[(ncb1 * somtab_dim2 + 1) *
6618 somtab_dim1], uvpav, isofav, ideriv, &somtab[(ncb1 *
6619 somtab_dim2 + 1) * somtab_dim1]);
6620 mma1noc_(uv11, &ibid1, &ii, &diftab[(ncb1 * diftab_dim2 + 1) *
6621 diftab_dim1], uvpav, isofav, ideriv, &diftab[(ncb1 *
6622 diftab_dim2 + 1) * diftab_dim1]);
6625 for (nd = 1; nd <= i__1; ++nd) {
6626 mma1noc_(uv11, &ncoeff[ncb1], &ii, &courbe[(nd + ncb1 *
6627 courbe_dim2) * courbe_dim1 + 1], uvpav, isofav, ideriv, &
6628 courbe[(nd + ncb1 * courbe_dim2) * courbe_dim1 + 1]);
6632 /* -> Update the nb of already created curves */
6635 /* -> ...otherwise try to cut the current interval in 2... */
6637 tmil = (tabdec[*nbcrbe + 1] + tabdec[*nbcrbe]) / 2.;
6640 ilong = (nupil - *nbcrbe) << 3;
6641 AdvApp2Var_SysBase::mcrfill_(&ilong, &tabdec[ideb],&tabdec[ideb1]);
6642 tabdec[ideb] = tmil;
6646 /* ---------- Make approximation of the rest -----------
6651 /* --------------------- Return code of error -----------------
6653 /* --> Pb with dynamic allocation */
6657 /* --> Inputs incoherent. */
6662 /* -------------------------- Dynamic desallocation -------------------
6667 anAdvApp2Var_SysBase.mcrdelt_(&c__8, &ndwrk, wrkar, &iofwr, &ier);
6674 /* ------------------------------ The end -------------------------------
6679 AdvApp2Var_SysBase::maermsg_("MMA2FNC", iercod, 7L);
6682 AdvApp2Var_SysBase::mgsomsg_("MMA2FNC", 7L);
6687 //=======================================================================
6688 //function : mma2fx6_
6690 //=======================================================================
6691 int AdvApp2Var_ApproxF2var::mma2fx6_(integer *ncfmxu,
6708 /* System generated locals */
6709 integer epsfro_dim1, epsfro_offset, patcan_dim1, patcan_dim2, patcan_dim3,
6710 patcan_dim4, patcan_offset, errmax_dim1, errmax_dim2,
6711 errmax_offset, ncoefu_dim1, ncoefu_offset, ncoefv_dim1,
6712 ncoefv_offset, i__1, i__2, i__3, i__4, i__5;
6713 doublereal d__1, d__2;
6715 /* Local variables */
6716 integer idim, ncfu, ncfv, id, ii, nd, jj, ku, kv, ns, ibb;
6720 /* **********************************************************************
6725 /* Reduction of degree when the squares are the squares of constraints. */
6729 /* TOUS,AB_SPECIFI::CARREAU&,REDUCTION,&CARREAU */
6731 /* INPUT ARGUMENTS : */
6732 /* ------------------ */
6733 /* NCFMXU: Max Nb of coeff by u of solution P(u,v) (table */
6734 /* PATCAN). This argument serves only to declare the size of this table. */
6735 /* NCFMXV: Max Nb of coeff by v of solution P(u,v) (table */
6736 /* PATCAN). This argument serves only to declare the size of this table. */
6737 /* NDIMEN: Total dimension of the space where the processed function */
6738 /* takes its values.(sum of dimensions of sub-spaces) */
6739 /* NBSESP: Nb of independent sub-spaces where the errors are measured. */
6740 /* NDIMSE: Table of dimensions of NBSESP sub-spaces. */
6741 /* NBUPAT: Nb of square solution by u. */
6742 /* NBVPAT: Nb of square solution by v. */
6743 /* IORDRU: Order of constraint imposed at the extremities of iso-V */
6744 /* = 0, the extremities of iso-V are calculated */
6745 /* = 1, additionally the 1st derivative in the direction of iso-V is calculated */
6746 /* = 2, additionally the 2nd derivative in the direction of iso-V is calculated */
6747 /* IORDRV: Ordre de contrainte impose aux extremites de l'iso-U */
6748 /* = 0, on calcule les extremites de l'iso-U. */
6749 /* = 1, additionally the 1st derivative in the direction of iso-U is calculated */
6750 /* = 2, additionally the 2nd derivative in the direction of iso-U is calculated */
6751 /* EPSAPR: Table of imposed precisions, sub-space by sub-space. */
6752 /* EPSFRO: Table of imposed precisions, sub-space by sub-space on the limits of squares. */
6753 /* PATCAN: Table of coeff. in the canonic base of squares P(u,v) calculated for (u,v) in (-1,1). */
6754 /* ERRMAX: Table of MAX errors (sub-space by sub-space) */
6755 /* committed in the approximation of F(u,v) by P(u,v). */
6756 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6757 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6759 /* OUTPUT ARGUMENTS : */
6760 /* ------------------- */
6761 /* NCOEFU: Table of Nb of significative coeffs. by u of calculated squares. */
6762 /* NCOEFV: Table of Nb of significative coeffs. by v of calculated squares. */
6764 /* COMMONS USED : */
6765 /* ---------------- */
6767 /* REFERENCES CALLED : */
6768 /* --------------------- */
6770 /* DESCRIPTION/NOTES/LIMITATIONS : */
6771 /* ------------------------------- */
6773 /* **********************************************************************
6776 /* Name of the routine */
6779 /* Parameter adjustments */
6780 epsfro_dim1 = *nbsesp;
6781 epsfro_offset = epsfro_dim1 * 5 + 1;
6782 epsfro -= epsfro_offset;
6785 ncoefv_dim1 = *nbupat;
6786 ncoefv_offset = ncoefv_dim1 + 1;
6787 ncoefv -= ncoefv_offset;
6788 ncoefu_dim1 = *nbupat;
6789 ncoefu_offset = ncoefu_dim1 + 1;
6790 ncoefu -= ncoefu_offset;
6791 errmax_dim1 = *nbsesp;
6792 errmax_dim2 = *nbupat;
6793 errmax_offset = errmax_dim1 * (errmax_dim2 + 1) + 1;
6794 errmax -= errmax_offset;
6795 patcan_dim1 = *ncfmxu;
6796 patcan_dim2 = *ncfmxv;
6797 patcan_dim3 = *ndimen;
6798 patcan_dim4 = *nbupat;
6799 patcan_offset = patcan_dim1 * (patcan_dim2 * (patcan_dim3 * (patcan_dim4
6801 patcan -= patcan_offset;
6804 ibb = AdvApp2Var_SysBase::mnfndeb_();
6806 AdvApp2Var_SysBase::mgenmsg_("MMA2FX6", 7L);
6811 for (jj = 1; jj <= i__1; ++jj) {
6813 for (ii = 1; ii <= i__2; ++ii) {
6814 ncfu = ncoefu[ii + jj * ncoefu_dim1];
6815 ncfv = ncoefv[ii + jj * ncoefv_dim1];
6817 /* ********************************************************************** */
6818 /* -------------------- Reduction of degree by U ------------------------- */
6819 /* ********************************************************************** */
6822 if (ncfu <= (*iordru + 1) << 1 && ncfu > 2) {
6826 for (ns = 1; ns <= i__3; ++ns) {
6829 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6830 tol = advapp_min(d__1,d__2);
6832 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6833 tol = advapp_min(d__1,d__2);
6835 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6836 tol = advapp_min(d__1,d__2);
6838 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6839 tol = advapp_min(d__1,d__2);
6840 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6843 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6844 tol = advapp_min(d__1,d__2);
6846 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6847 tol = advapp_min(d__1,d__2);
6849 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6850 tol = advapp_min(d__1,d__2);
6852 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6853 tol = advapp_min(d__1,d__2);
6858 for (nd = 1; nd <= i__4; ++nd) {
6861 for (kv = 1; kv <= i__5; ++kv) {
6862 bid += (d__1 = patcan[ncfu + (kv + (id + (ii + jj
6863 * patcan_dim4) * patcan_dim3) *
6864 patcan_dim2) * patcan_dim1], advapp_abs(d__1));
6870 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6871 errmax_dim2) * errmax_dim1]) {
6882 /* ********************************************************************** */
6883 /* -------------------- Reduction of degree by V ------------------------- */
6884 /* ********************************************************************** */
6887 if (ncfv <= (*iordrv + 1) << 1 && ncfv > 2) {
6891 for (ns = 1; ns <= i__3; ++ns) {
6894 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 9];
6895 tol = advapp_min(d__1,d__2);
6897 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 10];
6898 tol = advapp_min(d__1,d__2);
6900 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 11];
6901 tol = advapp_min(d__1,d__2);
6903 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 12];
6904 tol = advapp_min(d__1,d__2);
6905 if (ii == 1 || ii == *nbupat || jj == 1 || jj == *nbvpat)
6908 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 5];
6909 tol = advapp_min(d__1,d__2);
6911 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 6];
6912 tol = advapp_min(d__1,d__2);
6914 d__1 = tol, d__2 = epsfro[ns + epsfro_dim1 * 7];
6915 tol = advapp_min(d__1,d__2);
6917 d__1 = tol, d__2 = epsfro[ns + (epsfro_dim1 << 3)];
6918 tol = advapp_min(d__1,d__2);
6923 for (nd = 1; nd <= i__4; ++nd) {
6926 for (ku = 1; ku <= i__5; ++ku) {
6927 bid += (d__1 = patcan[ku + (ncfv + (id + (ii + jj
6928 * patcan_dim4) * patcan_dim3) *
6929 patcan_dim2) * patcan_dim1], advapp_abs(d__1));
6935 if (bid > tol * 1e-6 || bid > errmax[ns + (ii + jj *
6936 errmax_dim2) * errmax_dim1]) {
6947 /* --- Return the nbs of coeff. and pass to the next square --- */
6950 ncoefu[ii + jj * ncoefu_dim1] = advapp_max(ncfu,2);
6951 ncoefv[ii + jj * ncoefv_dim1] = advapp_max(ncfv,2);
6957 /* ------------------------------ The End -------------------------------
6961 AdvApp2Var_SysBase::mgsomsg_("MMA2FX6", 7L);
6967 //=======================================================================
6968 //function : mma2jmx_
6970 //=======================================================================
6971 int AdvApp2Var_ApproxF2var::mma2jmx_(integer *ndgjac,
6975 /* Initialized data */
6977 static doublereal xmax2[57] = { .9682458365518542212948163499456,
6978 .986013297183269340427888048593603,
6979 1.07810420343739860362585159028115,
6980 1.17325804490920057010925920756025,
6981 1.26476561266905634732910520370741,
6982 1.35169950227289626684434056681946,
6983 1.43424378958284137759129885012494,
6984 1.51281316274895465689402798226634,
6985 1.5878364329591908800533936587012,
6986 1.65970112228228167018443636171226,
6987 1.72874345388622461848433443013543,
6988 1.7952515611463877544077632304216,
6989 1.85947199025328260370244491818047,
6990 1.92161634324190018916351663207101,
6991 1.98186713586472025397859895825157,
6992 2.04038269834980146276967984252188,
6993 2.09730119173852573441223706382076,
6994 2.15274387655763462685970799663412,
6995 2.20681777186342079455059961912859,
6996 2.25961782459354604684402726624239,
6997 2.31122868752403808176824020121524,
6998 2.36172618435386566570998793688131,
6999 2.41117852396114589446497298177554,
7000 2.45964731268663657873849811095449,
7001 2.50718840313973523778244737914028,
7002 2.55385260994795361951813645784034,
7003 2.59968631659221867834697883938297,
7004 2.64473199258285846332860663371298,
7005 2.68902863641518586789566216064557,
7006 2.73261215675199397407027673053895,
7007 2.77551570192374483822124304745691,
7008 2.8177699459714315371037628127545,
7009 2.85940333797200948896046563785957,
7010 2.90044232019793636101516293333324,
7011 2.94091151970640874812265419871976,
7012 2.98083391718088702956696303389061,
7013 3.02023099621926980436221568258656,
7014 3.05912287574998661724731962377847,
7015 3.09752842783622025614245706196447,
7016 3.13546538278134559341444834866301,
7017 3.17295042316122606504398054547289,
7018 3.2099992681699613513775259670214,
7019 3.24662674946606137764916854570219,
7020 3.28284687953866689817670991319787,
7021 3.31867291347259485044591136879087,
7022 3.35411740487202127264475726990106,
7023 3.38919225660177218727305224515862,
7024 3.42390876691942143189170489271753,
7025 3.45827767149820230182596660024454,
7026 3.49230918177808483937957161007792,
7027 3.5260130200285724149540352829756,
7028 3.55939845146044235497103883695448,
7029 3.59247431368364585025958062194665,
7030 3.62524904377393592090180712976368,
7031 3.65773070318071087226169680450936,
7032 3.68992700068237648299565823810245,
7033 3.72184531357268220291630708234186 };
7034 static doublereal xmax4[55] = { 1.1092649593311780079813740546678,
7035 1.05299572648705464724876659688996,
7036 1.0949715351434178709281698645813,
7037 1.15078388379719068145021100764647,
7038 1.2094863084718701596278219811869,
7039 1.26806623151369531323304177532868,
7040 1.32549784426476978866302826176202,
7041 1.38142537365039019558329304432581,
7042 1.43575531950773585146867625840552,
7043 1.48850442653629641402403231015299,
7044 1.53973611681876234549146350844736,
7045 1.58953193485272191557448229046492,
7046 1.63797820416306624705258190017418,
7047 1.68515974143594899185621942934906,
7048 1.73115699602477936547107755854868,
7049 1.77604489805513552087086912113251,
7050 1.81989256661534438347398400420601,
7051 1.86276344480103110090865609776681,
7052 1.90471563564740808542244678597105,
7053 1.94580231994751044968731427898046,
7054 1.98607219357764450634552790950067,
7055 2.02556989246317857340333585562678,
7056 2.06433638992049685189059517340452,
7057 2.10240936014742726236706004607473,
7058 2.13982350649113222745523925190532,
7059 2.17661085564771614285379929798896,
7060 2.21280102016879766322589373557048,
7061 2.2484214321456956597803794333791,
7062 2.28349755104077956674135810027654,
7063 2.31805304852593774867640120860446,
7064 2.35210997297725685169643559615022,
7065 2.38568889602346315560143377261814,
7066 2.41880904328694215730192284109322,
7067 2.45148841120796359750021227795539,
7068 2.48374387161372199992570528025315,
7069 2.5155912654873773953959098501893,
7070 2.54704548720896557684101746505398,
7071 2.57812056037881628390134077704127,
7072 2.60882970619319538196517982945269,
7073 2.63918540521920497868347679257107,
7074 2.66919945330942891495458446613851,
7075 2.69888301230439621709803756505788,
7076 2.72824665609081486737132853370048,
7077 2.75730041251405791603760003778285,
7078 2.78605380158311346185098508516203,
7079 2.81451587035387403267676338931454,
7080 2.84269522483114290814009184272637,
7081 2.87060005919012917988363332454033,
7082 2.89823818258367657739520912946934,
7083 2.92561704377132528239806135133273,
7084 2.95274375377994262301217318010209,
7085 2.97962510678256471794289060402033,
7086 3.00626759936182712291041810228171,
7087 3.03267744830655121818899164295959,
7088 3.05886060707437081434964933864149 };
7089 static doublereal xmax6[53] = { 1.21091229812484768570102219548814,
7090 1.11626917091567929907256116528817,
7091 1.1327140810290884106278510474203,
7092 1.1679452722668028753522098022171,
7093 1.20910611986279066645602153641334,
7094 1.25228283758701572089625983127043,
7095 1.29591971597287895911380446311508,
7096 1.3393138157481884258308028584917,
7097 1.3821288728999671920677617491385,
7098 1.42420414683357356104823573391816,
7099 1.46546895108549501306970087318319,
7100 1.50590085198398789708599726315869,
7101 1.54550385142820987194251585145013,
7102 1.58429644271680300005206185490937,
7103 1.62230484071440103826322971668038,
7104 1.65955905239130512405565733793667,
7105 1.69609056468292429853775667485212,
7106 1.73193098017228915881592458573809,
7107 1.7671112206990325429863426635397,
7108 1.80166107681586964987277458875667,
7109 1.83560897003644959204940535551721,
7110 1.86898184653271388435058371983316,
7111 1.90180515174518670797686768515502,
7112 1.93410285411785808749237200054739,
7113 1.96589749778987993293150856865539,
7114 1.99721027139062501070081653790635,
7115 2.02806108474738744005306947877164,
7116 2.05846864831762572089033752595401,
7117 2.08845055210580131460156962214748,
7118 2.11802334209486194329576724042253,
7119 2.14720259305166593214642386780469,
7120 2.17600297710595096918495785742803,
7121 2.20443832785205516555772788192013,
7122 2.2325216999457379530416998244706,
7123 2.2602654243075083168599953074345,
7124 2.28768115912702794202525264301585,
7125 2.3147799369092684021274946755348,
7126 2.34157220782483457076721300512406,
7127 2.36806787963276257263034969490066,
7128 2.39427635443992520016789041085844,
7129 2.42020656255081863955040620243062,
7130 2.44586699364757383088888037359254,
7131 2.47126572552427660024678584642791,
7132 2.49641045058324178349347438430311,
7133 2.52130850028451113942299097584818,
7134 2.54596686772399937214920135190177,
7135 2.5703922285006754089328998222275,
7136 2.59459096001908861492582631591134,
7137 2.61856915936049852435394597597773,
7138 2.64233265984385295286445444361827,
7139 2.66588704638685848486056711408168,
7140 2.68923766976735295746679957665724,
7141 2.71238965987606292679677228666411 };
7143 /* System generated locals */
7146 /* Local variables */
7152 /* **********************************************************************
7157 /* Calculate the max of Jacobo polynoms multiplied by the weight on */
7158 /* (-1,1) for order 0,4,6 or Legendre. */
7162 /* LEGENDRE,APPROXIMATION,ERREUR. */
7164 /* INPUT ARGUMENTS : */
7165 /* ------------------ */
7166 /* NDGJAC: Nb of Jacobi coeff. of approximation. */
7167 /* IORDRE: Order of continuity (from -1 to 2) */
7169 /* OUTPUT ARGUMENTS : */
7170 /* ------------------- */
7171 /* XJACMX: Table of maximums of Jacobi polynoms. */
7173 /* COMMONS USED : */
7174 /* ---------------- */
7176 /* REFERENCES CALLED : */
7177 /* --------------------- */
7179 /* DESCRIPTION/NOTES/LIMITATIONS : */
7180 /* ----------------------------------- */
7183 /* ***********************************************************************
7185 /* Name of the routine */
7186 /* ----------------------------- Initialisations ------------------------
7189 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7191 AdvApp2Var_SysBase::mgenmsg_("MMA2JMX", 7L);
7194 numax = *ndgjac - ((*iordre + 1) << 1);
7195 if (*iordre == -1) {
7197 for (ii = 0; ii <= i__1; ++ii) {
7198 bid = (ii * 2. + 1.) / 2.;
7199 xjacmx[ii] = sqrt(bid);
7202 } else if (*iordre == 0) {
7204 for (ii = 0; ii <= i__1; ++ii) {
7205 xjacmx[ii] = xmax2[ii];
7208 } else if (*iordre == 1) {
7210 for (ii = 0; ii <= i__1; ++ii) {
7211 xjacmx[ii] = xmax4[ii];
7214 } else if (*iordre == 2) {
7216 for (ii = 0; ii <= i__1; ++ii) {
7217 xjacmx[ii] = xmax6[ii];
7222 /* ------------------------- The end ------------------------------------
7226 AdvApp2Var_SysBase::mgsomsg_("MMA2JMX", 7L);
7231 //=======================================================================
7232 //function : mma2moy_
7234 //=======================================================================
7235 int mma2moy_(integer *ndgumx,
7247 /* System generated locals */
7248 integer patjac_dim1, patjac_dim2, patjac_offset, i__1, i__2, i__3;
7250 /* Local variables */
7252 integer minu, minv, idebu, idebv, ii, nd, jj;
7253 doublereal bid0, bid1;
7256 /* **********************************************************************
7261 /* Calculate the average approximation error made when only */
7262 /* the coefficients of PATJAC of degree between */
7263 /* 2*(IORDRU+1) and MINDGU by U and 2*(IORDRV+1) and MINDGV by V are preserved. */
7267 /* LEGENDRE,APPROXIMATION, AVERAGE ERROR */
7269 /* INPUT ARGUMENTS : */
7270 /* ------------------ */
7271 /* NDGUMX: Dimension by U of table PATJAC. */
7272 /* NDGVMX: Dimension by V of table PATJAC. */
7273 /* NDIMEN: Dimension of the space. */
7274 /* MINDGU: Lower limit of the index by U of PATJAC coeff to be taken into account. */
7275 /* MAXDGU: Upper limit of the index by U of PATJAC coeff to be taken into account. */
7276 /* MINDGV: Lower limit of the index by V of PATJAC coeff to be taken into account. */
7277 /* MAXDGV: Upper limit of the index by V of PATJAC coeff to be taken into account. */
7278 /* IORDRU: Order of continuity by U provided by square PATJAC (from -1 to 2) */
7279 /* IORDRV: Order of continuity by V provided by square PATJAC (from -1 to 2) */
7280 /* PATJAC: Table of coeff. of the approximation square with */
7281 /* constraints of order IORDRU by U and IORDRV by V. */
7283 /* OUTPUT ARGUMENTS : */
7284 /* ------------------- */
7285 /* ERRMOY: Average error commited by preserving only the coeff of */
7286 /* PATJAC 2*(IORDRU+1) in MINDGU by U and 2*(IORDRV+1) in MINDGV by V. */
7288 /* COMMONS USED : */
7289 /* ---------------- */
7291 /* REFERENCES CALLED : */
7292 /* --------------------- */
7294 /* DESCRIPTION/NOTES/LIMITATIONS : */
7295 /* ----------------------------------- */
7296 /* Table PATJAC stores the coeff. Cij of */
7297 /* approximation square F(U,V). Indexes i and j show the degree by */
7298 /* U and by V of the base polynoms. These base polynoms are in the form: */
7300 /* ((1 - U*U)**(IORDRU+1)).J(i-2*(IORDRU+1)(U), where */
7302 /* polynom J(i-2*(IORDU+1)(U) is the Jacobi polynom of order */
7303 /* IORDRU+1 (the same by V by replacing U by V in the above expression). */
7305 /* The contribution to the average error of term Cij when */
7306 /* it is removed from PATJAC is Cij*Cij. */
7309 /* ***********************************************************************
7311 /* Name of the routine */
7314 /* ----------------------------- Initialisations ------------------------
7317 /* Parameter adjustments */
7318 patjac_dim1 = *ndgumx + 1;
7319 patjac_dim2 = *ndgvmx + 1;
7320 patjac_offset = patjac_dim1 * patjac_dim2;
7321 patjac -= patjac_offset;
7324 ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 3;
7326 AdvApp2Var_SysBase::mgenmsg_("MMA2MOY", 7L);
7329 idebu = (*iordru + 1) << 1;
7330 idebv = (*iordrv + 1) << 1;
7331 minu = advapp_max(idebu,*mindgu);
7332 minv = advapp_max(idebv,*mindgv);
7336 /* ------------------ Calculation of the upper bound of the average error ------------ */
7337 /* -------------------- when the coeff. of indexes from MINDGU to MAXDGU ------ */
7338 /* ---------------- by U and of indexes from MINDGV to MAXDGV by V are removed -------------- */
7341 for (nd = 1; nd <= i__1; ++nd) {
7343 for (jj = minv; jj <= i__2; ++jj) {
7345 for (ii = idebu; ii <= i__3; ++ii) {
7346 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7347 bid0 += bid1 * bid1;
7356 for (nd = 1; nd <= i__1; ++nd) {
7358 for (jj = idebv; jj <= i__2; ++jj) {
7360 for (ii = minu; ii <= i__3; ++ii) {
7361 bid1 = patjac[ii + (jj + nd * patjac_dim2) * patjac_dim1];
7362 bid0 += bid1 * bid1;
7370 /* ----------------------- Calculation of the average error -------------
7374 *errmoy = sqrt(bid0);
7376 /* ------------------------- The end ------------------------------------
7380 AdvApp2Var_SysBase::mgsomsg_("MMA2MOY", 7L);
7385 //=======================================================================
7386 //function : mma2roo_
7388 //=======================================================================
7389 int AdvApp2Var_ApproxF2var::mma2roo_(integer *nbpntu,
7394 /* System generated locals */
7397 /* Local variables */
7400 /* **********************************************************************
7405 /* Return roots of Legendre for discretisations. */
7409 /* TOUS, AB_SPECIFI::CONTRAINTE&, DISCRETISATION, &POINT */
7411 /* INPUT ARGUMENTS : */
7412 /* ------------------ */
7413 /* NBPNTU: Nb of INTERNAL parameters of discretization BY U. */
7414 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7415 /* NBPNTV: Nb of INTERNAL parameters of discretization BY V. */
7416 /* This is also the nb of root of the Legendre polynom where the discretization is done. */
7418 /* OUTPUT ARGUMENTS : */
7419 /* ------------------- */
7420 /* UROOTL: Table of parameters of discretisation ON (-1,1) BY U.
7422 /* VROOTL: Table of parameters of discretisation ON (-1,1) BY V.
7425 /* COMMONS USED : */
7426 /* ---------------- */
7428 /* REFERENCES CALLED : */
7429 /* --------------------- */
7431 /* DESCRIPTION/NOTES/LIMITATIONS : */
7432 /* ----------------------------------- */
7435 /* **********************************************************************
7438 /* Name of the routine */
7441 /* Parameter adjustments */
7446 ibb = AdvApp2Var_SysBase::mnfndeb_();
7448 AdvApp2Var_SysBase::mgenmsg_("MMA2ROO", 7L);
7451 /* ---------------- Return the POSITIVE roots on U ------------------
7454 AdvApp2Var_MathBase::mmrtptt_(nbpntu, &urootl[(*nbpntu + 1) / 2 + 1]);
7456 for (ii = 1; ii <= i__1; ++ii) {
7457 urootl[ii] = -urootl[*nbpntu - ii + 1];
7460 if (*nbpntu % 2 == 1) {
7461 urootl[*nbpntu / 2 + 1] = 0.;
7464 /* ---------------- Return the POSITIVE roots on V ------------------
7467 AdvApp2Var_MathBase::mmrtptt_(nbpntv, &vrootl[(*nbpntv + 1) / 2 + 1]);
7469 for (ii = 1; ii <= i__1; ++ii) {
7470 vrootl[ii] = -vrootl[*nbpntv - ii + 1];
7473 if (*nbpntv % 2 == 1) {
7474 vrootl[*nbpntv / 2 + 1] = 0.;
7477 /* ------------------------------ The End -------------------------------
7481 AdvApp2Var_SysBase::mgsomsg_("MMA2ROO", 7L);
7485 //=======================================================================
7486 //function : mmmapcoe_
7488 //=======================================================================
7489 int mmmapcoe_(integer *ndim,
7499 /* System generated locals */
7500 integer somtab_dim1, somtab_offset, diftab_dim1, diftab_offset,
7501 crvjac_dim1, crvjac_offset, gsstab_dim1, i__1, i__2, i__3;
7503 /* Local variables */
7504 integer igss, ikdeb;
7506 integer nd, ik, ir, nbroot, ibb;
7508 /* **********************************************************************
7513 /* Calculate the coefficients of polinomial approximation curve */
7514 /* of degree NDGJAC by the method of smallest squares starting from */
7515 /* the discretization of function on the roots of Legendre polynom */
7516 /* of degree NBPNTS. */
7520 /* FONCTION,APPROXIMATION,COEFFICIENT,POLYNOME */
7522 /* INPUT ARGUMENTS : */
7523 /* ------------------ */
7524 /* NDIM : Dimension of the space. */
7525 /* NDGJAC : Max Degree of the polynom of approximation. */
7526 /* The representation in the orthogonal base starts from degree */
7527 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7528 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7529 /* IORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7530 /* to step of constraints, C0,C1 or C2. */
7531 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7532 /* are calculated the coefficients of integration by */
7533 /* Gauss method. It is required to set NBPNTS=30,40,50 or 61 */
7534 /* and NDGJAC < NBPNTS. */
7535 /* SOMTAB : Table of F(ti)+F(-ti) with ti in ROOTAB. */
7536 /* DIFTAB : Table of F(ti)-F(-ti) with ti in ROOTAB. */
7537 /* GSSTAB(i,k) : Table of coefficients of integration by the Gauss method : */
7538 /* i varies from 0 to NBPNTS and */
7539 /* k varies from 0 to NDGJAC-2*(JORDRE+1). */
7541 /* OUTPUT ARGUMENTSE : */
7542 /* ------------------- */
7543 /* CRVJAC : Curve of approximation of FONCNP with eventually */
7544 /* taking into account of constraints at the extremities. */
7545 /* This curve is of degree NDGJAC. */
7547 /* COMMONS USED : */
7548 /* ---------------- */
7550 /* REFERENCES CALLED : */
7551 /* --------------------- */
7553 /* DESCRIPTION/NOTES/LIMITATIONS : */
7554 /* ------------------------------- */
7556 /* **********************************************************************
7559 /* Name of the routine */
7561 /* Parameter adjustments */
7562 crvjac_dim1 = *ndgjac + 1;
7563 crvjac_offset = crvjac_dim1;
7564 crvjac -= crvjac_offset;
7565 gsstab_dim1 = *nbpnts / 2 + 1;
7566 diftab_dim1 = *nbpnts / 2 + 1;
7567 diftab_offset = diftab_dim1;
7568 diftab -= diftab_offset;
7569 somtab_dim1 = *nbpnts / 2 + 1;
7570 somtab_offset = somtab_dim1;
7571 somtab -= somtab_offset;
7574 ibb = AdvApp2Var_SysBase::mnfndeb_();
7576 AdvApp2Var_SysBase::mgenmsg_("MMMAPCO", 7L);
7578 ikdeb = (*iordre + 1) << 1;
7579 nbroot = *nbpnts / 2;
7582 for (nd = 1; nd <= i__1; ++nd) {
7584 /* ----------------- Calculate the coefficients of even degree ----------
7588 for (ik = ikdeb; ik <= i__2; ik += 2) {
7592 for (ir = 1; ir <= i__3; ++ir) {
7593 bidon += somtab[ir + nd * somtab_dim1] * gsstab[ir + igss *
7597 crvjac[ik + nd * crvjac_dim1] = bidon;
7601 /* --------------- Calculate the coefficients of uneven degree ----------
7605 for (ik = ikdeb + 1; ik <= i__2; ik += 2) {
7609 for (ir = 1; ir <= i__3; ++ir) {
7610 bidon += diftab[ir + nd * diftab_dim1] * gsstab[ir + igss *
7614 crvjac[ik + nd * crvjac_dim1] = bidon;
7621 /* ------- Add terms connected to the supplementary root (0.D0) ------ */
7622 /* ----------- of Legendre polynom of uneven degree NBPNTS -----------
7625 if (*nbpnts % 2 == 0) {
7629 for (nd = 1; nd <= i__1; ++nd) {
7631 for (ik = ikdeb; ik <= i__2; ik += 2) {
7633 crvjac[ik + nd * crvjac_dim1] += somtab[nd * somtab_dim1] *
7634 gsstab[igss * gsstab_dim1];
7640 /* ------------------------------ The end -------------------------------
7645 AdvApp2Var_SysBase::mgsomsg_("MMMAPCO", 7L);
7649 //=======================================================================
7650 //function : mmaperm_
7652 //=======================================================================
7653 int mmaperm_(integer *ncofmx,
7661 /* System generated locals */
7662 integer crvjac_dim1, crvjac_offset, i__1, i__2;
7664 /* Local variables */
7666 integer i__, ia, nd, ncfcut, ibb;
7669 /* **********************************************************************
7674 /* Calculate the square root of the average quadratic error */
7675 /* of approximation done when only the */
7676 /* first NCFNEW coefficients of a curve of degree NCOEFF-1 */
7677 /* written in NORMALIZED Jacobi base of order 2*(IORDRE+1) are preserved. */
7681 /* LEGENDRE,POLYGONE,APPROXIMATION,ERREUR. */
7683 /* INPUT ARGUMENTS : */
7684 /* ------------------ */
7685 /* NCOFMX : Maximum degree of the curve. */
7686 /* NDIM : Dimension of the space. */
7687 /* NCOEFF : Degree +1 of the curve. */
7688 /* IORDRE : Order of constraint of continuity at the extremities. */
7689 /* CRVJAC : The curve the degree which of will be lowered. */
7690 /* NCFNEW : Degree +1 of the resulting polynom. */
7692 /* OUTPUT ARGUMENTS : */
7693 /* ------------------- */
7694 /* ERRMOY : Average precision of approximation. */
7696 /* COMMONS USED : */
7697 /* ---------------- */
7699 /* REFERENCES CALLED : */
7700 /* ----------------------- */
7702 /* DESCRIPTION/NOTES/LIMITATIONS : */
7703 /* ----------------------------------- */
7705 /* ***********************************************************************
7708 /* Name of the routine */
7710 /* Parameter adjustments */
7711 crvjac_dim1 = *ncofmx;
7712 crvjac_offset = crvjac_dim1 + 1;
7713 crvjac -= crvjac_offset;
7716 ibb = AdvApp2Var_SysBase::mnfndeb_();
7718 AdvApp2Var_SysBase::mgenmsg_("MMAPERM", 7L);
7721 /* --------- Minimum degree that can be reached : Stop at 1 or IA -------
7724 ia = (*iordre + 1) << 1;
7726 if (*ncfnew + 1 > ncfcut) {
7727 ncfcut = *ncfnew + 1;
7730 /* -------------- Elimination of coefficients of high degree ------------ */
7731 /* ----------- Loop on the series of Jacobi :NCFCUT --> NCOEFF --------- */
7736 for (nd = 1; nd <= i__1; ++nd) {
7738 for (i__ = ncfcut; i__ <= i__2; ++i__) {
7739 bidj = crvjac[i__ + nd * crvjac_dim1];
7746 /* ----------- Square Root of average quadratic error e -----------
7750 *errmoy = sqrt(bid);
7752 /* ------------------------------- The end ------------------------------
7756 AdvApp2Var_SysBase::mgsomsg_("MMAPERM", 7L);
7760 //=======================================================================
7761 //function : mmapptt_
7763 //=======================================================================
7764 int AdvApp2Var_ApproxF2var::mmapptt_(const integer *ndgjac,
7765 const integer *nbpnts,
7766 const integer *jordre,
7770 /* System generated locals */
7771 integer cgauss_dim1, i__1;
7773 /* Local variables */
7774 integer kjac, iptt, ipdb0, infdg, iptdb, mxjac, ilong, ibb;
7776 /* **********************************************************************
7781 /* Load the elements required for integration by */
7782 /* Gauss method to obtain the coefficients in the base of */
7783 /* Legendre of the approximation by the least squares of a */
7784 /* function. The elements are stored in commons MMAPGSS */
7785 /* (case without constraint), MMAPGS0 (constraints C0), MMAPGS1 */
7786 /* (constraints C1) and MMAPGS2 (constraints C2). */
7790 /* INTEGRATION,GAUSS,JACOBI */
7792 /* INPUT ARGUMENTS : */
7793 /* ------------------ */
7794 /* NDGJAC : Max degree of the polynom of approximation. */
7795 /* The representation in orthogonal base goes from degree */
7796 /* 0 to degree NDGJAC-2*(JORDRE+1). The polynomial base */
7797 /* is the base of Jacobi of order -1 (Legendre), 0, 1 and 2 */
7798 /* NBPNTS : Degree of the polynom of Legendre on the roots which of */
7799 /* are calculated the coefficients of integration by the */
7800 /* method of Gauss. It is required that NBPNTS=8,10,15,20,25, */
7801 /* 30,40,50 or 61 and NDGJAC < NBPNTS. */
7802 /* JORDRE : Order of the base of Jacobi (-1,0,1 or 2). Corresponds */
7803 /* to step of constraints C0,C1 or C2. */
7805 /* OUTPUT ARGUMENTS : */
7806 /* ------------------- */
7807 /* CGAUSS(i,k) : Table of coefficients of integration by */
7808 /* Gauss method : i varies from 0 to the integer part */
7809 /* of NBPNTS/2 and k varies from 0 to NDGJAC-2*(JORDRE+1). */
7810 /* These are the coeff. of integration associated to */
7811 /* positive roots of the polynom of Legendre of degree */
7812 /* NBPNTS. CGAUSS(0,k) contains coeff. */
7813 /* of integration associated to root t = 0 when */
7814 /* NBPNTS is uneven. */
7815 /* IERCOD : Error code. */
7817 /* = 11 NBPNTS is not 8,10,15,20,25,30,40,50 or 61. */
7818 /* = 21 JORDRE is not -1,0,1 or 2. */
7819 /* = 31 NDGJAC is too great or too small. */
7821 /* COMMONS USED : */
7822 /* ---------------- */
7823 /* MMAPGSS,MMAPGS0,MMAPGS1,MMAPGS2. */
7824 /* ***********************************************************************
7826 /* Parameter adjustments */
7827 cgauss_dim1 = *nbpnts / 2 + 1;
7830 ibb = AdvApp2Var_SysBase::mnfndeb_();
7832 AdvApp2Var_SysBase::mgenmsg_("MMAPPTT", 7L);
7836 /* ------------------- Tests on the validity of inputs ----------------
7839 infdg = (*jordre + 1) << 1;
7840 if (*nbpnts != 8 && *nbpnts != 10 && *nbpnts != 15 && *nbpnts != 20 && *
7841 nbpnts != 25 && *nbpnts != 30 && *nbpnts != 40 && *nbpnts != 50 &&
7846 if (*jordre < -1 || *jordre > 2) {
7850 if (*ndgjac >= *nbpnts || *ndgjac < infdg) {
7854 /* --------------- Calculation of the start pointer following NBPNTS -----------
7859 iptdb += (8 - infdg) << 2;
7862 iptdb += (10 - infdg) * 5;
7865 iptdb += (15 - infdg) * 7;
7868 iptdb += (20 - infdg) * 10;
7871 iptdb += (25 - infdg) * 12;
7874 iptdb += (30 - infdg) * 15;
7877 iptdb += (40 - infdg) * 20;
7880 iptdb += (50 - infdg) * 25;
7885 ipdb0 = ipdb0 + (14 - infdg) / 2 + 1;
7888 ipdb0 = ipdb0 + (24 - infdg) / 2 + 1;
7891 /* ------------------ Choice of the common depending on JORDRE -------------
7894 if (*jordre == -1) {
7907 /* ---------------- Common MMAPGSS (case without constraints) ----------------
7911 ilong = *nbpnts / 2 << 3;
7913 for (kjac = 0; kjac <= i__1; ++kjac) {
7914 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7915 AdvApp2Var_SysBase::mcrfill_(&ilong,
7916 &mmapgss_.gslxjs[iptt - 1],
7917 &cgauss[kjac * cgauss_dim1 + 1]);
7920 /* --> Case when the number of points is uneven. */
7921 if (*nbpnts % 2 == 1) {
7924 for (kjac = 0; kjac <= i__1; kjac += 2) {
7925 cgauss[kjac * cgauss_dim1] = mmapgss_.gsl0js[iptt - 1];
7930 for (kjac = 1; kjac <= i__1; kjac += 2) {
7931 cgauss[kjac * cgauss_dim1] = 0.;
7937 /* ---------------- Common MMAPGS0 (case with constraints C0) -------------
7941 mxjac = *ndgjac - infdg;
7942 ilong = *nbpnts / 2 << 3;
7944 for (kjac = 0; kjac <= i__1; ++kjac) {
7945 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7946 AdvApp2Var_SysBase::mcrfill_(&ilong,
7947 &mmapgs0_.gslxj0[iptt - 1],
7948 &cgauss[kjac * cgauss_dim1 + 1]);
7951 /* --> Case when the number of points is uneven. */
7952 if (*nbpnts % 2 == 1) {
7955 for (kjac = 0; kjac <= i__1; kjac += 2) {
7956 cgauss[kjac * cgauss_dim1] = mmapgs0_.gsl0j0[iptt - 1];
7961 for (kjac = 1; kjac <= i__1; kjac += 2) {
7962 cgauss[kjac * cgauss_dim1] = 0.;
7968 /* ---------------- Common MMAPGS1 (case with constraints C1) -------------
7972 mxjac = *ndgjac - infdg;
7973 ilong = *nbpnts / 2 << 3;
7975 for (kjac = 0; kjac <= i__1; ++kjac) {
7976 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
7977 AdvApp2Var_SysBase::mcrfill_(&ilong,
7978 &mmapgs1_.gslxj1[iptt - 1],
7979 &cgauss[kjac * cgauss_dim1 + 1]);
7982 /* --> Case when the number of points is uneven. */
7983 if (*nbpnts % 2 == 1) {
7986 for (kjac = 0; kjac <= i__1; kjac += 2) {
7987 cgauss[kjac * cgauss_dim1] = mmapgs1_.gsl0j1[iptt - 1];
7992 for (kjac = 1; kjac <= i__1; kjac += 2) {
7993 cgauss[kjac * cgauss_dim1] = 0.;
7999 /* ---------------- Common MMAPGS2 (case with constraints C2) -------------
8003 mxjac = *ndgjac - infdg;
8004 ilong = *nbpnts / 2 << 3;
8006 for (kjac = 0; kjac <= i__1; ++kjac) {
8007 iptt = iptdb + kjac * (*nbpnts / 2) + 1;
8008 AdvApp2Var_SysBase::mcrfill_(&ilong,
8009 &mmapgs2_.gslxj2[iptt - 1],
8010 &cgauss[kjac * cgauss_dim1 + 1]);
8013 /* --> Cas of uneven number of points. */
8014 if (*nbpnts % 2 == 1) {
8017 for (kjac = 0; kjac <= i__1; kjac += 2) {
8018 cgauss[kjac * cgauss_dim1] = mmapgs2_.gsl0j2[iptt - 1];
8023 for (kjac = 1; kjac <= i__1; kjac += 2) {
8024 cgauss[kjac * cgauss_dim1] = 0.;
8030 /* ------------------------- Return the error code --------------
8032 /* --> NBPNTS is not OK */
8036 /* --> JORDRE is not OK */
8040 /* --> NDGJAC is not OK */
8045 /* -------------------------------- The end -----------------------------
8050 AdvApp2Var_SysBase::maermsg_("MMAPPTT", iercod, 7L);
8053 AdvApp2Var_SysBase::mgsomsg_("MMAPPTT", 7L);
8059 //=======================================================================
8060 //function : mmjacpt_
8062 //=======================================================================
8063 int mmjacpt_(const integer *ndimen,
8064 const integer *ncoefu,
8065 const integer *ncoefv,
8066 const integer *iordru,
8067 const integer *iordrv,
8068 const doublereal *ptclgd,
8072 /* System generated locals */
8073 integer ptccan_dim1, ptccan_dim2, ptccan_offset, ptclgd_dim1, ptclgd_dim2,
8074 ptclgd_offset, ptcaux_dim1, ptcaux_dim2, ptcaux_dim3,
8075 ptcaux_offset, i__1, i__2, i__3;
8077 /* Local variables */
8078 integer kdim, nd, ii, jj, ibb;
8080 /* ***********************************************************************
8085 /* Passage from canonical to Jacobi base for a */
8086 /* "square" in a space of arbitrary dimension. */
8090 /* SMOOTHING,BASE,LEGENDRE */
8093 /* INPUT ARGUMENTS : */
8094 /* ------------------ */
8095 /* NDIMEN : Dimension of the space. */
8096 /* NCOEFU : Degree+1 by U. */
8097 /* NCOEFV : Degree+1 by V. */
8098 /* IORDRU : Order of Jacobi polynoms by U. */
8099 /* IORDRV : Order of Jacobi polynoms by V. */
8100 /* PTCLGD : The square in the Jacobi base. */
8102 /* OUTPUT ARGUMENTS : */
8103 /* ------------------- */
8104 /* PTCAUX : Auxilliary space. */
8105 /* PTCCAN : The square in the canonic base (-1,1) */
8107 /* COMMONS USED : */
8108 /* ---------------- */
8110 /* APPLIED REFERENCES : */
8111 /* ----------------------- */
8113 /* DESCRIPTION/NOTES/LIMITATIONS : */
8114 /* ----------------------------------- */
8115 /* Cancels and replaces MJACPC */
8117 /* *********************************************************************
8119 /* Name of the routine */
8122 /* Parameter adjustments */
8123 ptccan_dim1 = *ncoefu;
8124 ptccan_dim2 = *ncoefv;
8125 ptccan_offset = ptccan_dim1 * (ptccan_dim2 + 1) + 1;
8126 ptccan -= ptccan_offset;
8127 ptcaux_dim1 = *ncoefv;
8128 ptcaux_dim2 = *ncoefu;
8129 ptcaux_dim3 = *ndimen;
8130 ptcaux_offset = ptcaux_dim1 * (ptcaux_dim2 * (ptcaux_dim3 + 1) + 1) + 1;
8131 ptcaux -= ptcaux_offset;
8132 ptclgd_dim1 = *ncoefu;
8133 ptclgd_dim2 = *ncoefv;
8134 ptclgd_offset = ptclgd_dim1 * (ptclgd_dim2 + 1) + 1;
8135 ptclgd -= ptclgd_offset;
8138 ibb = AdvApp2Var_SysBase::mnfndeb_();
8140 AdvApp2Var_SysBase::mgenmsg_("MMJACPT", 7L);
8143 /* Passage into canonical by u. */
8145 kdim = *ndimen * *ncoefv;
8146 AdvApp2Var_MathBase::mmjaccv_(ncoefu,
8149 &ptclgd[ptclgd_offset],
8150 &ptcaux[ptcaux_offset],
8151 &ptccan[ptccan_offset]);
8153 /* Swapping of u and v. */
8156 for (nd = 1; nd <= i__1; ++nd) {
8158 for (jj = 1; jj <= i__2; ++jj) {
8160 for (ii = 1; ii <= i__3; ++ii) {
8161 ptcaux[jj + (ii + (nd + ptcaux_dim3) * ptcaux_dim2) *
8162 ptcaux_dim1] = ptccan[ii + (jj + nd * ptccan_dim2) *
8171 /* Passage into canonical by v. */
8173 kdim = *ndimen * *ncoefu;
8174 AdvApp2Var_MathBase::mmjaccv_(ncoefv,
8177 &ptcaux[((ptcaux_dim3 + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1],
8178 &ptccan[ptccan_offset],
8179 &ptcaux[(((ptcaux_dim3 << 1) + 1) * ptcaux_dim2 + 1) * ptcaux_dim1 + 1]);
8181 /* Swapping of u and v. */
8184 for (nd = 1; nd <= i__1; ++nd) {
8186 for (jj = 1; jj <= i__2; ++jj) {
8188 for (ii = 1; ii <= i__3; ++ii) {
8189 ptccan[ii + (jj + nd * ptccan_dim2) * ptccan_dim1] = ptcaux[
8190 jj + (ii + (nd + (ptcaux_dim3 << 1)) * ptcaux_dim2) *
8199 /* ---------------------------- THAT'S ALL FOLKS ------------------------
8203 AdvApp2Var_SysBase::mgsomsg_("MMJACPT", 7L);