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 } |