Doxygen Source Code Documentation
cl1.c File Reference
#include <math.h>#include <stddef.h>#include <stdio.h>#include <stdlib.h>#include "f2c.h"Go to the source code of this file.
Functions | |
| int | cl1_fort (integer *k, integer *l, integer *m, integer *n, integer *klmd, integer *klm2d, integer *nklmd, integer *n2d, real *q, integer *kode, real *toler, integer *iter, real *x, real *res, real *error, real *cu, integer *iu, integer *s) |
| float | cl1_solve (int ndim, int nvec, float *z, float **A, float *y, int cony) |
| float | cl1_solve_res (int ndim, int nvec, float *z, float **A, float *y, int cony, float *rez, int conr) |
Function Documentation
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 389 of file cl1.c. References abs, dabs, l, n1, n2, and q. Referenced by cl1_solve(), and cl1_solve_res().
00393 {
00394 /* System generated locals */
00395 integer q_dim1, q_offset, i__1, i__2;
00396 real r__1;
00397
00398 /* Local variables */
00399 static integer iimn, nklm;
00400 static real xmin, xmax;
00401 static integer iout, i__, j;
00402 static real z__;
00403 static integer iineg, maxit, n1, n2;
00404 static real pivot;
00405 static integer ia, ii, kk, in, nk, js;
00406 static real sn;
00407 static integer iphase, kforce;
00408 static real zu, zv;
00409 static integer nk1;
00410 static real tpivot;
00411 static integer klm, jmn, nkl, jpn;
00412 static real cuv;
00413 static doublereal sum;
00414 static integer klm1, klm2, nkl1;
00415
00416
00417 /* THIS SUBROUTINE USES A MODIFICATION OF THE SIMPLEX */
00418 /* METHOD OF LINEAR PROGRAMMING TO CALCULATE AN L1 SOLUTION */
00419 /* TO A K BY N SYSTEM OF LINEAR EQUATIONS */
00420 /* AX=B */
00421 /* SUBJECT TO L LINEAR EQUALITY CONSTRAINTS */
00422 /* CX=D */
00423 /* AND M LINEAR INEQUALITY CONSTRAINTS */
00424 /* EX.LE.F. */
00425
00426 /* DESCRIPTION OF PARAMETERS */
00427
00428 /* K NUMBER OF ROWS OF THE MATRIX A (K.GE.1). */
00429 /* L NUMBER OF ROWS OF THE MATRIX C (L.GE.0). */
00430 /* M NUMBER OF ROWS OF THE MATRIX E (M.GE.0). */
00431 /* N NUMBER OF COLUMNS OF THE MATRICES A,C,E (N.GE.1). */
00432 /* KLMD SET TO AT LEAST K+L+M FOR ADJUSTABLE DIMENSIONS. */
00433 /* KLM2D SET TO AT LEAST K+L+M+2 FOR ADJUSTABLE DIMENSIONS. */
00434 /* NKLMD SET TO AT LEAST N+K+L+M FOR ADJUSTABLE DIMENSIONS. */
00435 /* N2D SET TO AT LEAST N+2 FOR ADJUSTABLE DIMENSIONS */
00436
00437 /* Q TWO DIMENSIONAL REAL ARRAY WITH KLM2D ROWS AND */
00438 /* AT LEAST N2D COLUMNS. */
00439 /* ON ENTRY THE MATRICES A,C AND E, AND THE VECTORS */
00440 /* B,D AND F MUST BE STORED IN THE FIRST K+L+M ROWS */
00441 /* AND N+1 COLUMNS OF Q AS FOLLOWS */
00442 /* A B */
00443 /* Q = C D */
00444 /* E F */
00445 /* THESE VALUES ARE DESTROYED BY THE SUBROUTINE. */
00446
00447 /* KODE A CODE USED ON ENTRY TO, AND EXIT */
00448 /* FROM, THE SUBROUTINE. */
00449 /* ON ENTRY, THIS SHOULD NORMALLY BE SET TO 0. */
00450 /* HOWEVER, IF CERTAIN NONNEGATIVITY CONSTRAINTS */
00451 /* ARE TO BE INCLUDED IMPLICITLY, RATHER THAN */
00452 /* EXPLICITLY IN THE CONSTRAINTS EX.LE.F, THEN KODE */
00453 /* SHOULD BE SET TO 1, AND THE NONNEGATIVITY */
00454 /* CONSTRAINTS INCLUDED IN THE ARRAYS X AND */
00455 /* RES (SEE BELOW). */
00456 /* ON EXIT, KODE HAS ONE OF THE */
00457 /* FOLLOWING VALUES */
00458 /* 0- OPTIMAL SOLUTION FOUND, */
00459 /* 1- NO FEASIBLE SOLUTION TO THE */
00460 /* CONSTRAINTS, */
00461 /* 2- CALCULATIONS TERMINATED */
00462 /* PREMATURELY DUE TO ROUNDING ERRORS, */
00463 /* 3- MAXIMUM NUMBER OF ITERATIONS REACHED. */
00464
00465 /* TOLER A SMALL POSITIVE TOLERANCE. EMPIRICAL */
00466 /* EVIDENCE SUGGESTS TOLER = 10**(-D*2/3), */
00467 /* WHERE D REPRESENTS THE NUMBER OF DECIMAL */
00468 /* DIGITS OF ACCURACY AVAILABLE. ESSENTIALLY, */
00469 /* THE SUBROUTINE CANNOT DISTINGUISH BETWEEN ZERO */
00470 /* AND ANY QUANTITY WHOSE MAGNITUDE DOES NOT EXCEED */
00471 /* TOLER. IN PARTICULAR, IT WILL NOT PIVOT ON ANY */
00472 /* NUMBER WHOSE MAGNITUDE DOES NOT EXCEED TOLER. */
00473
00474 /* ITER ON ENTRY ITER MUST CONTAIN AN UPPER BOUND ON */
00475 /* THE MAXIMUM NUMBER OF ITERATIONS ALLOWED. */
00476 /* A SUGGESTED VALUE IS 10*(K+L+M). ON EXIT ITER */
00477 /* GIVES THE NUMBER OF SIMPLEX ITERATIONS. */
00478
00479 /* X ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST N2D. */
00480 /* ON EXIT THIS ARRAY CONTAINS A */
00481 /* SOLUTION TO THE L1 PROBLEM. IF KODE=1 */
00482 /* ON ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE */
00483 /* SIMPLE NONNEGATIVITY CONSTRAINTS ON THE */
00484 /* VARIABLES. THE VALUES -1, 0, OR 1 */
00485 /* FOR X(J) INDICATE THAT THE J-TH VARIABLE */
00486 /* IS RESTRICTED TO BE .LE.0, UNRESTRICTED, */
00487 /* OR .GE.0 RESPECTIVELY. */
00488
00489 /* RES ONE DIMENSIONAL REAL ARRAY OF SIZE AT LEAST KLMD. */
00490 /* ON EXIT THIS CONTAINS THE RESIDUALS B-AX */
00491 /* IN THE FIRST K COMPONENTS, D-CX IN THE */
00492 /* NEXT L COMPONENTS (THESE WILL BE =0),AND */
00493 /* F-EX IN THE NEXT M COMPONENTS. IF KODE=1 ON */
00494 /* ENTRY, THIS ARRAY IS ALSO USED TO INCLUDE SIMPLE */
00495 /* NONNEGATIVITY CONSTRAINTS ON THE RESIDUALS */
00496 /* B-AX. THE VALUES -1, 0, OR 1 FOR RES(I) */
00497 /* INDICATE THAT THE I-TH RESIDUAL (1.LE.I.LE.K) IS */
00498 /* RESTRICTED TO BE .LE.0, UNRESTRICTED, OR .GE.0 */
00499 /* RESPECTIVELY. */
00500
00501 /* ERROR ON EXIT, THIS GIVES THE MINIMUM SUM OF */
00502 /* ABSOLUTE VALUES OF THE RESIDUALS. */
00503
00504 /* CU A TWO DIMENSIONAL REAL ARRAY WITH TWO ROWS AND */
00505 /* AT LEAST NKLMD COLUMNS USED FOR WORKSPACE. */
00506
00507 /* IU A TWO DIMENSIONAL INTEGER ARRAY WITH TWO ROWS AND */
00508 /* AT LEAST NKLMD COLUMNS USED FOR WORKSPACE. */
00509
00510 /* S INTEGER ARRAY OF SIZE AT LEAST KLMD, USED FOR */
00511 /* WORKSPACE. */
00512
00513 /* IF YOUR FORTRAN COMPILER PERMITS A SINGLE COLUMN OF A TWO */
00514 /* DIMENSIONAL ARRAY TO BE PASSED TO A ONE DIMENSIONAL ARRAY */
00515 /* THROUGH A SUBROUTINE CALL, CONSIDERABLE SAVINGS IN */
00516 /* EXECUTION TIME MAY BE ACHIEVED THROUGH THE USE OF THE */
00517 /* FOLLOWING SUBROUTINE, WHICH OPERATES ON COLUMN VECTORS. */
00518 /* SUBROUTINE COL(V1, V2, XMLT, NOTROW, K) */
00519 /* THIS SUBROUTINE ADDS TO THE VECTOR V1 A MULTIPLE OF THE */
00520 /* VECTOR V2 (ELEMENTS 1 THROUGH K EXCLUDING NOTROW). */
00521 /* DIMENSION V1(K), V2(K) */
00522 /* KEND = NOTROW - 1 */
00523 /* KSTART = NOTROW + 1 */
00524 /* IF (KEND .LT. 1) GO TO 20 */
00525 /* DO 10 I=1,KEND */
00526 /* V1(I) = V1(I) + XMLT*V2(I) */
00527 /* 10 CONTINUE */
00528 /* IF(KSTART .GT. K) GO TO 40 */
00529 /* 20 DO 30 I=KSTART,K */
00530 /* V1(I) = V1(I) + XMLT*V2(I) */
00531 /* 30 CONTINUE */
00532 /* 40 RETURN */
00533 /* END */
00534 /* SEE COMMENTS FOLLOWING STATEMENT LABELLED 440 FOR */
00535 /* INSTRUCTIONS ON THE IMPLEMENTATION OF THIS MODIFICATION. */
00536 /* -----------------------------------------------------------------------
00537 */
00538
00539 /* INITIALIZATION. */
00540
00541 /* Parameter adjustments */
00542 --s;
00543 --res;
00544 iu -= 3;
00545 cu -= 3;
00546 --x;
00547 q_dim1 = *klm2d;
00548 q_offset = q_dim1 + 1;
00549 q -= q_offset;
00550
00551 /* Function Body */
00552 maxit = *iter;
00553 n1 = *n + 1;
00554 n2 = *n + 2;
00555 nk = *n + *k;
00556 nk1 = nk + 1;
00557 nkl = nk + *l;
00558 nkl1 = nkl + 1;
00559 klm = *k + *l + *m;
00560 klm1 = klm + 1;
00561 klm2 = klm + 2;
00562 nklm = *n + klm;
00563 kforce = 1;
00564 *iter = 0;
00565 js = 1;
00566 ia = 0;
00567 /* SET UP LABELS IN Q. */
00568 i__1 = *n;
00569 for (j = 1; j <= i__1; ++j) {
00570 q[klm2 + j * q_dim1] = (real) j;
00571 /* L10: */
00572 }
00573 i__1 = klm;
00574 for (i__ = 1; i__ <= i__1; ++i__) {
00575 q[i__ + n2 * q_dim1] = (real) (*n + i__);
00576 if (q[i__ + n1 * q_dim1] >= 0.f) {
00577 goto L30;
00578 }
00579 i__2 = n2;
00580 for (j = 1; j <= i__2; ++j) {
00581 q[i__ + j * q_dim1] = -q[i__ + j * q_dim1];
00582 /* L20: */
00583 }
00584 L30:
00585 ;
00586 }
00587 /* SET UP PHASE 1 COSTS. */
00588 iphase = 2;
00589 i__1 = nklm;
00590 for (j = 1; j <= i__1; ++j) {
00591 cu[(j << 1) + 1] = 0.f;
00592 cu[(j << 1) + 2] = 0.f;
00593 iu[(j << 1) + 1] = 0;
00594 iu[(j << 1) + 2] = 0;
00595 /* L40: */
00596 }
00597 if (*l == 0) {
00598 goto L60;
00599 }
00600 i__1 = nkl;
00601 for (j = nk1; j <= i__1; ++j) {
00602 cu[(j << 1) + 1] = 1.f;
00603 cu[(j << 1) + 2] = 1.f;
00604 iu[(j << 1) + 1] = 1;
00605 iu[(j << 1) + 2] = 1;
00606 /* L50: */
00607 }
00608 iphase = 1;
00609 L60:
00610 if (*m == 0) {
00611 goto L80;
00612 }
00613 i__1 = nklm;
00614 for (j = nkl1; j <= i__1; ++j) {
00615 cu[(j << 1) + 2] = 1.f;
00616 iu[(j << 1) + 2] = 1;
00617 jmn = j - *n;
00618 if (q[jmn + n2 * q_dim1] < 0.f) {
00619 iphase = 1;
00620 }
00621 /* L70: */
00622 }
00623 L80:
00624 if (*kode == 0) {
00625 goto L150;
00626 }
00627 i__1 = *n;
00628 for (j = 1; j <= i__1; ++j) {
00629 if ((r__1 = x[j]) < 0.f) {
00630 goto L90;
00631 } else if (r__1 == 0) {
00632 goto L110;
00633 } else {
00634 goto L100;
00635 }
00636 L90:
00637 cu[(j << 1) + 1] = 1.f;
00638 iu[(j << 1) + 1] = 1;
00639 goto L110;
00640 L100:
00641 cu[(j << 1) + 2] = 1.f;
00642 iu[(j << 1) + 2] = 1;
00643 L110:
00644 ;
00645 }
00646 i__1 = *k;
00647 for (j = 1; j <= i__1; ++j) {
00648 jpn = j + *n;
00649 if ((r__1 = res[j]) < 0.f) {
00650 goto L120;
00651 } else if (r__1 == 0) {
00652 goto L140;
00653 } else {
00654 goto L130;
00655 }
00656 L120:
00657 cu[(jpn << 1) + 1] = 1.f;
00658 iu[(jpn << 1) + 1] = 1;
00659 if (q[j + n2 * q_dim1] > 0.f) {
00660 iphase = 1;
00661 }
00662 goto L140;
00663 L130:
00664 cu[(jpn << 1) + 2] = 1.f;
00665 iu[(jpn << 1) + 2] = 1;
00666 if (q[j + n2 * q_dim1] < 0.f) {
00667 iphase = 1;
00668 }
00669 L140:
00670 ;
00671 }
00672 L150:
00673 if (iphase == 2) {
00674 goto L500;
00675 }
00676 /* COMPUTE THE MARGINAL COSTS. */
00677 L160:
00678 i__1 = n1;
00679 for (j = js; j <= i__1; ++j) {
00680 sum = 0.;
00681 i__2 = klm;
00682 for (i__ = 1; i__ <= i__2; ++i__) {
00683 ii = q[i__ + n2 * q_dim1];
00684 if (ii < 0) {
00685 goto L170;
00686 }
00687 z__ = cu[(ii << 1) + 1];
00688 goto L180;
00689 L170:
00690 iineg = -ii;
00691 z__ = cu[(iineg << 1) + 2];
00692 L180:
00693 sum += (doublereal) q[i__ + j * q_dim1] * (doublereal) z__;
00694 /* L190: */
00695 }
00696 q[klm1 + j * q_dim1] = sum;
00697 /* L200: */
00698 }
00699 i__1 = *n;
00700 for (j = js; j <= i__1; ++j) {
00701 ii = q[klm2 + j * q_dim1];
00702 if (ii < 0) {
00703 goto L210;
00704 }
00705 z__ = cu[(ii << 1) + 1];
00706 goto L220;
00707 L210:
00708 iineg = -ii;
00709 z__ = cu[(iineg << 1) + 2];
00710 L220:
00711 q[klm1 + j * q_dim1] -= z__;
00712 /* L230: */
00713 }
00714 /* DETERMINE THE VECTOR TO ENTER THE BASIS. */
00715 L240:
00716 xmax = 0.f;
00717 if (js > *n) {
00718 goto L490;
00719 }
00720 i__1 = *n;
00721 for (j = js; j <= i__1; ++j) {
00722 zu = q[klm1 + j * q_dim1];
00723 ii = q[klm2 + j * q_dim1];
00724 if (ii > 0) {
00725 goto L250;
00726 }
00727 ii = -ii;
00728 zv = zu;
00729 zu = -zu - cu[(ii << 1) + 1] - cu[(ii << 1) + 2];
00730 goto L260;
00731 L250:
00732 zv = -zu - cu[(ii << 1) + 1] - cu[(ii << 1) + 2];
00733 L260:
00734 if (kforce == 1 && ii > *n) {
00735 goto L280;
00736 }
00737 if (iu[(ii << 1) + 1] == 1) {
00738 goto L270;
00739 }
00740 if (zu <= xmax) {
00741 goto L270;
00742 }
00743 xmax = zu;
00744 in = j;
00745 L270:
00746 if (iu[(ii << 1) + 2] == 1) {
00747 goto L280;
00748 }
00749 if (zv <= xmax) {
00750 goto L280;
00751 }
00752 xmax = zv;
00753 in = j;
00754 L280:
00755 ;
00756 }
00757 if (xmax <= *toler) {
00758 goto L490;
00759 }
00760 if (q[klm1 + in * q_dim1] == xmax) {
00761 goto L300;
00762 }
00763 i__1 = klm2;
00764 for (i__ = 1; i__ <= i__1; ++i__) {
00765 q[i__ + in * q_dim1] = -q[i__ + in * q_dim1];
00766 /* L290: */
00767 }
00768 q[klm1 + in * q_dim1] = xmax;
00769 /* DETERMINE THE VECTOR TO LEAVE THE BASIS. */
00770 L300:
00771 if (iphase == 1 || ia == 0) {
00772 goto L330;
00773 }
00774 xmax = 0.f;
00775 i__1 = ia;
00776 for (i__ = 1; i__ <= i__1; ++i__) {
00777 z__ = (r__1 = q[i__ + in * q_dim1], dabs(r__1));
00778 if (z__ <= xmax) {
00779 goto L310;
00780 }
00781 xmax = z__;
00782 iout = i__;
00783 L310:
00784 ;
00785 }
00786 if (xmax <= *toler) {
00787 goto L330;
00788 }
00789 i__1 = n2;
00790 for (j = 1; j <= i__1; ++j) {
00791 z__ = q[ia + j * q_dim1];
00792 q[ia + j * q_dim1] = q[iout + j * q_dim1];
00793 q[iout + j * q_dim1] = z__;
00794 /* L320: */
00795 }
00796 iout = ia;
00797 --ia;
00798 pivot = q[iout + in * q_dim1];
00799 goto L420;
00800 L330:
00801 kk = 0;
00802 i__1 = klm;
00803 for (i__ = 1; i__ <= i__1; ++i__) {
00804 z__ = q[i__ + in * q_dim1];
00805 if (z__ <= *toler) {
00806 goto L340;
00807 }
00808 ++kk;
00809 res[kk] = q[i__ + n1 * q_dim1] / z__;
00810 s[kk] = i__;
00811 L340:
00812 ;
00813 }
00814 L350:
00815 if (kk > 0) {
00816 goto L360;
00817 }
00818 *kode = 2;
00819 goto L590;
00820 L360:
00821 xmin = res[1];
00822 iout = s[1];
00823 j = 1;
00824 if (kk == 1) {
00825 goto L380;
00826 }
00827 i__1 = kk;
00828 for (i__ = 2; i__ <= i__1; ++i__) {
00829 if (res[i__] >= xmin) {
00830 goto L370;
00831 }
00832 j = i__;
00833 xmin = res[i__];
00834 iout = s[i__];
00835 L370:
00836 ;
00837 }
00838 res[j] = res[kk];
00839 s[j] = s[kk];
00840 L380:
00841 --kk;
00842 pivot = q[iout + in * q_dim1];
00843 ii = q[iout + n2 * q_dim1];
00844 if (iphase == 1) {
00845 goto L400;
00846 }
00847 if (ii < 0) {
00848 goto L390;
00849 }
00850 if (iu[(ii << 1) + 2] == 1) {
00851 goto L420;
00852 }
00853 goto L400;
00854 L390:
00855 iineg = -ii;
00856 if (iu[(iineg << 1) + 1] == 1) {
00857 goto L420;
00858 }
00859 L400:
00860 ii = abs(ii);
00861 cuv = cu[(ii << 1) + 1] + cu[(ii << 1) + 2];
00862 if (q[klm1 + in * q_dim1] - pivot * cuv <= *toler) {
00863 goto L420;
00864 }
00865 /* BYPASS INTERMEDIATE VERTICES. */
00866 i__1 = n1;
00867 for (j = js; j <= i__1; ++j) {
00868 z__ = q[iout + j * q_dim1];
00869 q[klm1 + j * q_dim1] -= z__ * cuv;
00870 q[iout + j * q_dim1] = -z__;
00871 /* L410: */
00872 }
00873 q[iout + n2 * q_dim1] = -q[iout + n2 * q_dim1];
00874 goto L350;
00875 /* GAUSS-JORDAN ELIMINATION. */
00876 L420:
00877 if (*iter < maxit) {
00878 goto L430;
00879 }
00880 *kode = 3;
00881 goto L590;
00882 L430:
00883 ++(*iter);
00884 i__1 = n1;
00885 for (j = js; j <= i__1; ++j) {
00886 if (j != in) {
00887 q[iout + j * q_dim1] /= pivot;
00888 }
00889 /* L440: */
00890 }
00891 /* IF PERMITTED, USE SUBROUTINE COL OF THE DESCRIPTION */
00892 /* SECTION AND REPLACE THE FOLLOWING SEVEN STATEMENTS DOWN */
00893 /* TO AND INCLUDING STATEMENT NUMBER 460 BY.. */
00894 /* DO 460 J=JS,N1 */
00895 /* IF(J .EQ. IN) GO TO 460 */
00896 /* Z = -Q(IOUT,J) */
00897 /* CALL COL(Q(1,J), Q(1,IN), Z, IOUT, KLM1) */
00898 /* 460 CONTINUE */
00899 i__1 = n1;
00900 for (j = js; j <= i__1; ++j) {
00901 if (j == in) {
00902 goto L460;
00903 }
00904 z__ = -q[iout + j * q_dim1];
00905 i__2 = klm1;
00906 for (i__ = 1; i__ <= i__2; ++i__) {
00907 if (i__ != iout) {
00908 q[i__ + j * q_dim1] += z__ * q[i__ + in * q_dim1];
00909 }
00910 /* L450: */
00911 }
00912 L460:
00913 ;
00914 }
00915 tpivot = -pivot;
00916 i__1 = klm1;
00917 for (i__ = 1; i__ <= i__1; ++i__) {
00918 if (i__ != iout) {
00919 q[i__ + in * q_dim1] /= tpivot;
00920 }
00921 /* L470: */
00922 }
00923 q[iout + in * q_dim1] = 1.f / pivot;
00924 z__ = q[iout + n2 * q_dim1];
00925 q[iout + n2 * q_dim1] = q[klm2 + in * q_dim1];
00926 q[klm2 + in * q_dim1] = z__;
00927 ii = dabs(z__);
00928 if (iu[(ii << 1) + 1] == 0 || iu[(ii << 1) + 2] == 0) {
00929 goto L240;
00930 }
00931 i__1 = klm2;
00932 for (i__ = 1; i__ <= i__1; ++i__) {
00933 z__ = q[i__ + in * q_dim1];
00934 q[i__ + in * q_dim1] = q[i__ + js * q_dim1];
00935 q[i__ + js * q_dim1] = z__;
00936 /* L480: */
00937 }
00938 ++js;
00939 goto L240;
00940 /* TEST FOR OPTIMALITY. */
00941 L490:
00942 if (kforce == 0) {
00943 goto L580;
00944 }
00945 if (iphase == 1 && q[klm1 + n1 * q_dim1] <= *toler) {
00946 goto L500;
00947 }
00948 kforce = 0;
00949 goto L240;
00950 /* SET UP PHASE 2 COSTS. */
00951 L500:
00952 iphase = 2;
00953 i__1 = nklm;
00954 for (j = 1; j <= i__1; ++j) {
00955 cu[(j << 1) + 1] = 0.f;
00956 cu[(j << 1) + 2] = 0.f;
00957 /* L510: */
00958 }
00959 i__1 = nk;
00960 for (j = n1; j <= i__1; ++j) {
00961 cu[(j << 1) + 1] = 1.f;
00962 cu[(j << 1) + 2] = 1.f;
00963 /* L520: */
00964 }
00965 i__1 = klm;
00966 for (i__ = 1; i__ <= i__1; ++i__) {
00967 ii = q[i__ + n2 * q_dim1];
00968 if (ii > 0) {
00969 goto L530;
00970 }
00971 ii = -ii;
00972 if (iu[(ii << 1) + 2] == 0) {
00973 goto L560;
00974 }
00975 cu[(ii << 1) + 2] = 0.f;
00976 goto L540;
00977 L530:
00978 if (iu[(ii << 1) + 1] == 0) {
00979 goto L560;
00980 }
00981 cu[(ii << 1) + 1] = 0.f;
00982 L540:
00983 ++ia;
00984 i__2 = n2;
00985 for (j = 1; j <= i__2; ++j) {
00986 z__ = q[ia + j * q_dim1];
00987 q[ia + j * q_dim1] = q[i__ + j * q_dim1];
00988 q[i__ + j * q_dim1] = z__;
00989 /* L550: */
00990 }
00991 L560:
00992 ;
00993 }
00994 goto L160;
00995 L570:
00996 if (q[klm1 + n1 * q_dim1] <= *toler) {
00997 goto L500;
00998 }
00999 *kode = 1;
01000 goto L590;
01001 L580:
01002 if (iphase == 1) {
01003 goto L570;
01004 }
01005 /* PREPARE OUTPUT. */
01006 *kode = 0;
01007 L590:
01008 sum = 0.;
01009 i__1 = *n;
01010 for (j = 1; j <= i__1; ++j) {
01011 x[j] = 0.f;
01012 /* L600: */
01013 }
01014 i__1 = klm;
01015 for (i__ = 1; i__ <= i__1; ++i__) {
01016 res[i__] = 0.f;
01017 /* L610: */
01018 }
01019 i__1 = klm;
01020 for (i__ = 1; i__ <= i__1; ++i__) {
01021 ii = q[i__ + n2 * q_dim1];
01022 sn = 1.f;
01023 if (ii > 0) {
01024 goto L620;
01025 }
01026 ii = -ii;
01027 sn = -1.f;
01028 L620:
01029 if (ii > *n) {
01030 goto L630;
01031 }
01032 x[ii] = sn * q[i__ + n1 * q_dim1];
01033 goto L640;
01034 L630:
01035 iimn = ii - *n;
01036 res[iimn] = sn * q[i__ + n1 * q_dim1];
01037 if (ii >= n1 && ii <= nk) {
01038 sum += (doublereal) q[i__ + n1 * q_dim1];
01039 }
01040 L640:
01041 ;
01042 }
01043 *error = sum;
01044 return 0;
01045 } /* cl1_ */
|
|
||||||||||||||||||||||||||||
|
Definition at line 37 of file cl1.c. References calloc, cl1_fort(), free, l, and q. Referenced by L1F_worker(), and main().
00038 {
00039 /* loop counters */
00040
00041 int jj , ii ;
00042
00043 /* variables for CL1 (types declared in f2c.h) */
00044
00045 integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode=0,iter , *iu,*s ;
00046 real *q , toler , *x , *res , error , *cu ;
00047
00048 /*-- check inputs --*/
00049
00050 if( ndim < 1 || nvec < 1 ) kode = 4 ;
00051 if( A == NULL || y == NULL || z == NULL ) kode = 4 ;
00052 for( jj=0 ; jj < nvec ; jj++ ) if( A[jj] == NULL ) kode = 4 ;
00053
00054 if( kode ){
00055 fprintf(stderr,"** cl1_solve ERROR: illegal inputs!\n") ;
00056 return (float)(-kode) ;
00057 }
00058
00059 /*-- setup call to CL1 --*/
00060
00061 k = ndim ;
00062 l = 0 ; /* no linear equality constraints */
00063 m = 0 ; /* no linear inequality constraints */
00064 n = nvec ;
00065
00066 klmd = k+l+m ;
00067 klm2d = k+l+m+2 ;
00068 nklmd = n+k+l+m ;
00069 n2d = n+2 ;
00070
00071 kode = (cony != 0) ; /* enforce implicit constraints on x[] */
00072 iter = 11*klmd ;
00073
00074 toler = 0.0001 ;
00075 error = 0.0 ;
00076
00077 /* input/output matrices & vectors */
00078
00079 q = (real *) calloc( 1, sizeof(real) * klm2d*n2d ) ;
00080 x = (real *) calloc( 1, sizeof(real) * n2d ) ;
00081 res = (real *) calloc( 1, sizeof(real) * klmd ) ;
00082
00083 /* workspaces */
00084
00085 cu = (real *) calloc( 1, sizeof(real) * 2*nklmd ) ;
00086 iu = (integer *) calloc( 1, sizeof(integer) * 2*nklmd ) ;
00087 s = (integer *) calloc( 1, sizeof(integer) * klmd ) ;
00088
00089 /* load matrices & vectors */
00090
00091 for( jj=0 ; jj < nvec ; jj++ )
00092 for( ii=0 ; ii < ndim ; ii++ )
00093 q[ii+jj*klm2d] = A[jj][ii] ; /* matrix */
00094
00095 for( ii=0 ; ii < ndim ; ii++ )
00096 q[ii+nvec*klm2d] = z[ii] ; /* vector */
00097
00098 if( cony ){
00099 for( jj=0 ; jj < nvec ; jj++ ) /* constraints on solution */
00100 x[jj] = (y[jj] < 0.0) ? -1.0
00101 :(y[jj] > 0.0) ? 1.0 : 0.0 ;
00102 }
00103
00104 for( ii=0 ; ii < klmd ; ii++ ) /* no constraints on resids */
00105 res[ii] = 0.0 ;
00106
00107 /*-- do the work --*/
00108
00109 cl1_fort( &k, &l, &m, &n,
00110 &klmd, &klm2d, &nklmd, &n2d,
00111 q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
00112
00113 free(q) ; free(res) ; free(cu) ; free(iu) ; free(s) ;
00114
00115 if( kode != 0 ){
00116 free(x) ;
00117 #if 0
00118 switch( kode ){
00119 case 1: fprintf(stderr,"** cl1_solve ERROR: no feasible solution!\n"); break;
00120 case 2: fprintf(stderr,"** cl1_solve ERROR: rounding errors!\n") ; break;
00121 case 3: fprintf(stderr,"** cl1_solve ERROR: max iterations!\n") ; break;
00122 default: fprintf(stderr,"** cl1_solve ERROR: unknown problem!\n") ; break;
00123 }
00124 #endif
00125 return (float)(-kode) ;
00126 }
00127
00128 /*-- copy results into output --*/
00129
00130 for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ;
00131
00132 free(x) ; return (float)error ;
00133 }
|
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 170 of file cl1.c. References calloc, cl1_fort(), free, l, and q. Referenced by main().
00172 {
00173 /* loop counters */
00174
00175 int jj , ii ;
00176
00177 /* variables for CL1 (types declared in f2c.h) */
00178
00179 integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode=0,iter , *iu,*s ;
00180 real *q , toler , *x , *res , error , *cu ;
00181
00182 /*-- check inputs --*/
00183
00184 if( ndim < 1 || nvec < 1 ) kode = 4 ;
00185 if( A == NULL || y == NULL || z == NULL ) kode = 4 ;
00186 for( jj=0 ; jj < nvec ; jj++ ) if( A[jj] == NULL ) kode = 4 ;
00187
00188 if( kode ){
00189 fprintf(stderr,"** cl1_solve ERROR: illegal inputs!\n") ;
00190 return (float)(-kode) ;
00191 }
00192
00193 /*-- setup call to CL1 --*/
00194
00195 k = ndim ;
00196 l = 0 ; /* no linear equality constraints */
00197 m = 0 ; /* no linear inequality constraints */
00198 n = nvec ;
00199
00200 klmd = k+l+m ;
00201 klm2d = k+l+m+2 ;
00202 nklmd = n+k+l+m ;
00203 n2d = n+2 ;
00204
00205 kode = (cony != 0 || conr != 0) ; /* enforce implicit constraints on x[] */
00206 iter = 11*klmd ;
00207
00208 toler = 0.0001 ;
00209 error = 0.0 ;
00210
00211 /* input/output matrices & vectors */
00212
00213 q = (real *) calloc( 1, sizeof(real) * klm2d*n2d ) ;
00214 x = (real *) calloc( 1, sizeof(real) * n2d ) ;
00215 res = (real *) calloc( 1, sizeof(real) * klmd ) ;
00216
00217 /* workspaces */
00218
00219 cu = (real *) calloc( 1, sizeof(real) * 2*nklmd ) ;
00220 iu = (integer *) calloc( 1, sizeof(integer) * 2*nklmd ) ;
00221 s = (integer *) calloc( 1, sizeof(integer) * klmd ) ;
00222
00223 /* load matrices & vectors */
00224
00225 for( jj=0 ; jj < nvec ; jj++ )
00226 for( ii=0 ; ii < ndim ; ii++ )
00227 q[ii+jj*klm2d] = A[jj][ii] ; /* matrix */
00228
00229 for( ii=0 ; ii < ndim ; ii++ )
00230 q[ii+nvec*klm2d] = z[ii] ; /* vector */
00231
00232 if( cony ){
00233 for( jj=0 ; jj < nvec ; jj++ ) /* constraints on solution */
00234 x[jj] = (y[jj] < 0.0) ? -1.0
00235 :(y[jj] > 0.0) ? 1.0 : 0.0 ;
00236 }
00237
00238 if( conr ){
00239 for( ii=0 ; ii < ndim ; ii++ ) /* constraints on resids */
00240 res[ii] = (rez[ii] < 0.0) ? -1.0
00241 :(rez[ii] > 0.0) ? 1.0 : 0.0 ;
00242 }
00243
00244 /*-- do the work --*/
00245
00246 cl1_fort( &k, &l, &m, &n,
00247 &klmd, &klm2d, &nklmd, &n2d,
00248 q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
00249
00250 free(q) ; free(cu) ; free(iu) ; free(s) ;
00251
00252 if( kode != 0 ){
00253 free(x) ; free(res) ;
00254 #if 0
00255 switch( kode ){
00256 case 1: fprintf(stderr,"** cl1_solve ERROR: no feasible solution!\n"); break;
00257 case 2: fprintf(stderr,"** cl1_solve ERROR: rounding errors!\n") ; break;
00258 case 3: fprintf(stderr,"** cl1_solve ERROR: max iterations!\n") ; break;
00259 default: fprintf(stderr,"** cl1_solve ERROR: unknown problem!\n") ; break;
00260 }
00261 #endif
00262 return (float)(-kode) ;
00263 }
00264
00265 /*-- copy results into output --*/
00266
00267 for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ;
00268
00269 for( ii=0 ; ii < ndim ; ii++ ) rez[ii] = (float) res[ii] ;
00270
00271 free(res); free(x); return (float)error;
00272 }
|
