00001 #include <math.h>
00002 #include <stddef.h>
00003 #include <stdio.h>
00004 #include <stdlib.h>
00005 #include "f2c.h"
00006
00007
00008
00009 static int cl1_fort(integer *k, integer *l, integer *m, integer *n,
00010 integer *klmd, integer *klm2d, integer *nklmd,
00011 integer *n2d, real *q,
00012 integer *kode, real *toler, integer *iter, real *x,
00013 real *res, real * error, real *cu, integer *iu, integer *s) ;
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 float cl1_solve( int ndim, int nvec, float *z, float **A, float *y, int cony )
00038 {
00039
00040
00041 int jj , ii ;
00042
00043
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
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
00060
00061 k = ndim ;
00062 l = 0 ;
00063 m = 0 ;
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) ;
00072 iter = 11*klmd ;
00073
00074 toler = 0.0001 ;
00075 error = 0.0 ;
00076
00077
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
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
00090
00091 for( jj=0 ; jj < nvec ; jj++ )
00092 for( ii=0 ; ii < ndim ; ii++ )
00093 q[ii+jj*klm2d] = A[jj][ii] ;
00094
00095 for( ii=0 ; ii < ndim ; ii++ )
00096 q[ii+nvec*klm2d] = z[ii] ;
00097
00098 if( cony ){
00099 for( jj=0 ; jj < nvec ; jj++ )
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++ )
00105 res[ii] = 0.0 ;
00106
00107
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
00129
00130 for( jj=0 ; jj < nvec ; jj++ ) y[jj] = (float) x[jj] ;
00131
00132 free(x) ; return (float)error ;
00133 }
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170 float cl1_solve_res( int ndim, int nvec, float *z, float **A,
00171 float *y, int cony , float *rez , int conr )
00172 {
00173
00174
00175 int jj , ii ;
00176
00177
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
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
00194
00195 k = ndim ;
00196 l = 0 ;
00197 m = 0 ;
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) ;
00206 iter = 11*klmd ;
00207
00208 toler = 0.0001 ;
00209 error = 0.0 ;
00210
00211
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
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
00224
00225 for( jj=0 ; jj < nvec ; jj++ )
00226 for( ii=0 ; ii < ndim ; ii++ )
00227 q[ii+jj*klm2d] = A[jj][ii] ;
00228
00229 for( ii=0 ; ii < ndim ; ii++ )
00230 q[ii+nvec*klm2d] = z[ii] ;
00231
00232 if( cony ){
00233 for( jj=0 ; jj < nvec ; jj++ )
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++ )
00240 res[ii] = (rez[ii] < 0.0) ? -1.0
00241 :(rez[ii] > 0.0) ? 1.0 : 0.0 ;
00242 }
00243
00244
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
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 }
00273
00274 #if 0
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301 int cl1_pos_sum( int ndim , int nvec ,
00302 float * base_vec , float ** extra_vec , float * ec )
00303 {
00304
00305
00306 int jj , ii ;
00307
00308
00309
00310 integer k,l,m,n,klmd,klm2d,nklmd,n2d , kode,iter , *iu,*s ;
00311 real *q , toler , *x , *res , error , *cu ;
00312
00313
00314
00315 if( ndim < 1 || nvec < 1 ) return 4 ;
00316 if( base_vec == NULL || extra_vec == NULL || ec == NULL ) return 4 ;
00317 for( jj=0 ; jj < nvec ; jj++ ) if( extra_vec[jj] == NULL ) return 4 ;
00318
00319
00320
00321 k = ndim ;
00322 l = 0 ;
00323 m = 0 ;
00324 n = nvec ;
00325
00326 klmd = k+l+m ;
00327 klm2d = k+l+m+2 ;
00328 nklmd = n+k+l+m ;
00329 n2d = n+2 ;
00330
00331 kode = 1 ;
00332 iter = 10*klmd ;
00333
00334 toler = 0.0001 ;
00335 error = 0.0 ;
00336
00337
00338
00339 q = (real *) malloc( sizeof(real) * klm2d*n2d ) ;
00340 x = (real *) malloc( sizeof(real) * n2d ) ;
00341 res = (real *) malloc( sizeof(real) * klmd ) ;
00342
00343
00344
00345 cu = (real *) malloc( sizeof(real) * 2*nklmd ) ;
00346 iu = (integer *) malloc( sizeof(integer) * 2*nklmd ) ;
00347 s = (integer *) malloc( sizeof(integer) * klmd ) ;
00348
00349
00350
00351 for( jj=0 ; jj < nvec ; jj++ )
00352 for( ii=0 ; ii < ndim ; ii++ )
00353 q[ii+jj*klm2d] = extra_vec[jj][ii] ;
00354
00355 for( ii=0 ; ii < ndim ; ii++ )
00356 q[ii+nvec*klm2d] = base_vec[ii] ;
00357
00358 for( jj=0 ; jj < nvec ; jj++ ) x[jj] = 0.0 ;
00359
00360 for( ii=0 ; ii < ndim ; ii++ ) res[ii] = 1.0 ;
00361
00362
00363
00364 cl1_fort( &k, &l, &m, &n,
00365 &klmd, &klm2d, &nklmd, &n2d,
00366 q, &kode, &toler, &iter, x, res, &error, cu, iu, s) ;
00367
00368 free(q) ; free(res) ; free(cu) ; free(iu) ; free(s) ;
00369
00370 if( kode != 0 ){ free(x) ; return (int)kode ; }
00371
00372
00373
00374 for( jj=0 ; jj < nvec ; jj++ ) ec[jj] = (float)(-x[jj]) ;
00375
00376 free(x) ; return 0 ;
00377 }
00378 #endif
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389 static int cl1_fort(integer *k, integer *l, integer *m, integer *n,
00390 integer *klmd, integer *klm2d, integer *nklmd, integer *n2d, real *q,
00391 integer *kode, real *toler, integer *iter, real *x, real *res, real *
00392 error, real *cu, integer *iu, integer *s)
00393 {
00394
00395 integer q_dim1, q_offset, i__1, i__2;
00396 real r__1;
00397
00398
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
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
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
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
00568 i__1 = *n;
00569 for (j = 1; j <= i__1; ++j) {
00570 q[klm2 + j * q_dim1] = (real) j;
00571
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
00583 }
00584 L30:
00585 ;
00586 }
00587
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
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
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
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
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
00695 }
00696 q[klm1 + j * q_dim1] = sum;
00697
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
00713 }
00714
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
00767 }
00768 q[klm1 + in * q_dim1] = xmax;
00769
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
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
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
00872 }
00873 q[iout + n2 * q_dim1] = -q[iout + n2 * q_dim1];
00874 goto L350;
00875
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
00890 }
00891
00892
00893
00894
00895
00896
00897
00898
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
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
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
00937 }
00938 ++js;
00939 goto L240;
00940
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
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
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
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
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
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
01013 }
01014 i__1 = klm;
01015 for (i__ = 1; i__ <= i__1; ++i__) {
01016 res[i__] = 0.f;
01017
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 }
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471