luksan/mssubs.c

   1 #include <math.h>
   2 #include "luksan.h"
   3
   4 #define max(a,b) ((a) > (b) ? (a) : (b))
   5 #define iabs(a) ((a) < 0 ? -(a) : (a))
   6
   7 /*     subroutines extracted from mssubs.for */
   8 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
   9 /* FUNCTION MXVMAX             ALL SYSTEMS                   91/12/01 */
  10 /* PURPOSE : */
  11 /* L-INFINITY NORM OF A VECTOR. */
  12
  13 /* PARAMETERS : */
  14 /*  II  N  VECTOR DIMENSION. */
  15 /*  RI  X(N)  INPUT VECTOR. */
  16 /*  RR  MXVMAX  L-INFINITY NORM OF THE VECTOR X. */
  17
  18 double luksan_mxvmax__(int *n, double *x)
  19 {
  20     /* System generated locals */
  21     int i__1;
  22     double d__1, d__2, d__3;
  23
  24     /* Local variables */
  25     int i__;
  26     double mxvmax;
  27
  28     /* Parameter adjustments */
  29     --x;
  30
  31     /* Function Body */
  32     mxvmax = 0.;
  33     i__1 = *n;
  34     for (i__ = 1; i__ <= i__1; ++i__) {
  35 /* Computing MAX */
  36         d__2 = mxvmax, d__3 = (d__1 = x[i__], fabs(d__1));
  37         mxvmax = max(d__2,d__3);
  38 /* L1: */
  39     }
  40     return mxvmax;
  41 } /* luksan_mxvmax__ */
  42
  43 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
  44 /* SUBROUTINE MXVINE             ALL SYSTEMS                   94/12/01 */
  45 /* PURPOSE : */
  46 /* ELEMENTS OF THE INTEGER VECTOR ARE REPLACED BY THEIR ABSOLUTE VALUES. */
  47
  48 /* PARAMETERS : */
  49 /*  II  N DIMENSION OF THE INTEGER VECTOR. */
  50 /*  IU  IX(N)  INTEGER VECTOR WHICH IS UPDATED SO THAT IX(I):=ABS(IX(I)) */
  51 /*         FOR ALL I. */
  52
  53 void luksan_mxvine__(int *n, int *ix)
  54 {
  55     /* System generated locals */
  56     int i__1, i__2;
  57
  58     /* Local variables */
  59     int i__;
  60
  61     /* Parameter adjustments */
  62     --ix;
  63
  64     /* Function Body */
  65     i__1 = *n;
  66     for (i__ = 1; i__ <= i__1; ++i__) {
  67         ix[i__] = (i__2 = ix[i__], iabs(i__2));
  68 /* L1: */
  69     }
  70     return;
  71 } /* luksan_mxvine__ */
  72
  73 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
  74 /* SUBROUTINE MXDCMV               ALL SYSTEMS                91/12/01 */
  75 /* PURPOSE : */
  76 /* RANK-TWO UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A. */
  77 /* THIS MATRIX IS UPDATED BY THE RULE A:=A+ALF*X*TRANS(U)+BET*Y*TRANS(V). */
  78
  79 /* PARAMETERS : */
  80 /*  II  N  NUMBER OF ROWS OF THE MATRIX A. */
  81 /*  II  M  NUMBER OF COLUMNS OF THE MATRIX A. */
  82 /*  RU  A(N*M)  RECTANGULAR MATRIX STORED COLUMNWISE IN THE */
  83 /*         ONE-DIMENSIONAL ARRAY. */
  84 /*  RI  ALF  SCALAR PARAMETER. */
  85 /*  RI  X(N)  INPUT VECTOR. */
  86 /*  RI  U(M)  INPUT VECTOR. */
  87 /*  RI  BET  SCALAR PARAMETER. */
  88 /*  RI  Y(N)  INPUT VECTOR. */
  89 /*  RI  V(M)  INPUT VECTOR. */
  90
  91 void luksan_mxdcmv__(int *n, int *m, double *a,
  92         double *alf, double *x, double *u, double *bet,
  93         double *y, double *v)
  94 {
  95     /* System generated locals */
  96     int i__1, i__2;
  97
  98     /* Local variables */
  99     int i__, j, k;
 100     double tempa, tempb;
 101
 102     /* Parameter adjustments */
 103     --v;
 104     --y;
 105     --u;
 106     --x;
 107     --a;
 108
 109     /* Function Body */
 110     k = 0;
 111     i__1 = *m;
 112     for (j = 1; j <= i__1; ++j) {
 113         tempa = *alf * u[j];
 114         tempb = *bet * v[j];
 115         i__2 = *n;
 116         for (i__ = 1; i__ <= i__2; ++i__) {
 117             a[k + i__] = a[k + i__] + tempa * x[i__] + tempb * y[i__];
 118 /* L1: */
 119         }
 120         k += *n;
 121 /* L2: */
 122     }
 123     return;
 124 } /* luksan_mxdcmv__ */
 125
 126 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 127 /* SUBROUTINE MXVSAV                ALL SYSTEMS                91/12/01 */
 128 /* PORTABILITY : ALL SYSTEMS */
 129 /* 91/12/01 LU : ORIGINAL VERSION */
 130
 131 /* PURPOSE : */
 132 /* DIFFERENCE OF TWO VECTORS RETURNED IN THE SUBSTRACTED ONE. */
 133
 134 /* PARAMETERS : */
 135 /*  II  N  VECTOR DIMENSION. */
 136 /*  RI  X(N)  INPUT VECTOR. */
 137 /*  RU  Y(N)  UPDATE VECTOR WHERE Y:= X - Y. */
 138
 139 void luksan_mxvsav__(int *n, double *x, double *y)
 140 {
 141     /* System generated locals */
 142     int i__1;
 143
 144     /* Local variables */
 145     int i__;
 146     double temp;
 147
 148     /* Parameter adjustments */
 149     --y;
 150     --x;
 151
 152     /* Function Body */
 153     i__1 = *n;
 154     for (i__ = 1; i__ <= i__1; ++i__) {
 155         temp = y[i__];
 156         y[i__] = x[i__] - y[i__];
 157         x[i__] = temp;
 158 /* L10: */
 159     }
 160     return;
 161 } /* luksan_mxvsav__ */
 162
 163 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 164 /* SUBROUTINE MXVLIN                ALL SYSTEMS                92/12/01 */
 165 /* PURPOSE : */
 166 /* LINEAR COMBINATION OF TWO VECTORS. */
 167
 168 /* PARAMETERS : */
 169 /*  II  N  VECTOR DIMENSION. */
 170 /*  RI  A  SCALING FACTOR. */
 171 /*  RI  X(N)  INPUT VECTOR. */
 172 /*  RI  B  SCALING FACTOR. */
 173 /*  RI  Y(N)  INPUT VECTOR. */
 174 /*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= A*X + B*Y. */
 175
 176 void luksan_mxvlin__(int *n, double *a, double *x,
 177         double *b, double *y, double *z__)
 178 {
 179     /* System generated locals */
 180     int i__1;
 181
 182     /* Local variables */
 183     int i__;
 184
 185     /* Parameter adjustments */
 186     --z__;
 187     --y;
 188     --x;
 189
 190     /* Function Body */
 191     i__1 = *n;
 192     for (i__ = 1; i__ <= i__1; ++i__) {
 193         z__[i__] = *a * x[i__] + *b * y[i__];
 194 /* L1: */
 195     }
 196     return;
 197 } /* luksan_mxvlin__ */
 198
 199 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 200 /* SUBROUTINE MXDCMU               ALL SYSTEMS                91/12/01 */
 201 /* PURPOSE : */
 202 /* UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A. THIS MATRIX */
 203 /* IS UPDATED BY THE RULE A:=A+ALF*X*TRANS(Y). */
 204
 205 /* PARAMETERS : */
 206 /*  II  N  NUMBER OF ROWS OF THE MATRIX A. */
 207 /*  II  M  NUMBER OF COLUMNS OF THE MATRIX A. */
 208 /*  RU  A(N*M)  RECTANGULAR MATRIX STORED COLUMNWISE IN THE */
 209 /*         ONE-DIMENSIONAL ARRAY. */
 210 /*  RI  ALF  SCALAR PARAMETER. */
 211 /*  RI  X(N)  INPUT VECTOR. */
 212 /*  RI  Y(M)  INPUT VECTOR. */
 213
 214 void luksan_mxdcmu__(int *n, int *m, double *a,
 215         double *alf, double *x, double *y)
 216 {
 217     /* System generated locals */
 218     int i__1, i__2;
 219
 220     /* Local variables */
 221     int i__, j, k;
 222     double temp;
 223
 224     /* Parameter adjustments */
 225     --y;
 226     --x;
 227     --a;
 228
 229     /* Function Body */
 230     k = 0;
 231     i__1 = *m;
 232     for (j = 1; j <= i__1; ++j) {
 233         temp = *alf * y[j];
 234         i__2 = *n;
 235         for (i__ = 1; i__ <= i__2; ++i__) {
 236             a[k + i__] += temp * x[i__];
 237 /* L1: */
 238         }
 239         k += *n;
 240 /* L2: */
 241     }
 242     return;
 243 } /* luksan_mxdcmu__ */
 244
 245 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 246 /* SUBROUTINE MXVDIR                ALL SYSTEMS                91/12/01 */
 247 /* PURPOSE : */
 248 /* VECTOR AUGMENTED BY THE SCALED VECTOR. */
 249
 250 /* PARAMETERS : */
 251 /*  II  N  VECTOR DIMENSION. */
 252 /*  RI  A  SCALING FACTOR. */
 253 /*  RI  X(N)  INPUT VECTOR. */
 254 /*  RI  Y(N)  INPUT VECTOR. */
 255 /*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= Y + A*X. */
 256
 257 void luksan_mxvdir__(int *n, double *a, double *x,
 258                     double *y, double *z__)
 259 {
 260     /* System generated locals */
 261     int i__1;
 262
 263     /* Local variables */
 264     int i__;
 265
 266     /* Parameter adjustments */
 267     --z__;
 268     --y;
 269     --x;
 270
 271     /* Function Body */
 272     i__1 = *n;
 273     for (i__ = 1; i__ <= i__1; ++i__) {
 274         z__[i__] = y[i__] + *a * x[i__];
 275 /* L10: */
 276     }
 277 } /* luksan_mxvdir__ */
 278
 279 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 280 /* SUBROUTINE MXDCMD               ALL SYSTEMS                91/12/01
 281 * PURPOSE :
 282 * MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A
 283 * BY A VECTOR X AND ADDITION OF THE SCALED VECTOR ALF*Y.
 284 *
 285 * PARAMETERS :
 286 *  II  N  NUMBER OF ROWS OF THE MATRIX A.
 287 *  II  M  NUMBER OF COLUMNS OF THE MATRIX A.
 288 *  RI  A(N*M)  RECTANGULAR MATRIX STORED COLUMNWISE IN THE
 289 *         ONE-DIMENSIONAL ARRAY.
 290 *  RI  X(M)  INPUT VECTOR.
 291 *  RI  ALF  SCALING FACTOR.
 292 *  RI  Y(N)  INPUT VECTOR.
 293 *  RO  Z(N)  OUTPUT VECTOR EQUAL TO A*X+ALF*Y.
 294 *
 295 * SUBPROGRAMS USED :
 296 *  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
 297 *  S   MXVSCL  SCALING OF A VECTOR.
 298 */
 299 void luksan_mxdcmd__(int *n, int *m, double *a,
 300         double *x, double *alf, double *y, double *z__)
 301 {
 302     /* System generated locals */
 303     int i__1;
 304
 305     /* Local variables */
 306     int j, k;
 307
 308     /* Parameter adjustments */
 309     --z__;
 310     --y;
 311     --x;
 312     --a;
 313
 314     /* Function Body */
 315     luksan_mxvscl__(n, alf, &y[1], &z__[1]);
 316     k = 0;
 317     i__1 = *m;
 318     for (j = 1; j <= i__1; ++j) {
 319         luksan_mxvdir__(n, &x[j], &a[k + 1], &z__[1], &z__[1]);
 320         k += *n;
 321 /* L1: */
 322     }
 323     return;
 324 } /* luksan_mxdcmd__ */
 325
 326 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 327 /* SUBROUTINE MXDRCB               ALL SYSTEMS                91/12/01
 328 * PURPOSE :
 329 * BACKWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION OF
 330 * THE VECTOR X BY AN IMPLICIT BFGS UPDATE.
 331 *
 332 * PARAMETERS :
 333 *  II  N  NUMBER OF ROWS OF THE MATRICES A AND B.
 334 *  II  M  NUMBER OF COLUMNS OF THE MATRICES A AND B.
 335 *  RI  A(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 336 *  RI  B(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 337 *  RI  U(M)  VECTOR OF SCALAR COEFFICIENTS.
 338 *  RO  V(M)  VECTOR OF SCALAR COEFFICIENTS.
 339 *  RU  X(N)  PREMULTIPLIED VECTOR.
 340 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
 341 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
 342 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
 343 *         IX(I).EQ.-5.
 344 *
 345 * SUBPROGRAM USED :
 346 *  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
 347 *  RF  MXUDOT  DOT PRODUCT OF VECTORS.
 348 *
 349 * METHOD :
 350 * H.MATTHIES, G.STRANG: THE SOLUTION OF NONLINEAR FINITE ELEMENT
 351 * EQUATIONS. INT.J.NUMER. METHODS ENGN. 14 (1979) 1613-1626.
 352 */
 353 void luksan_mxdrcb__(int *n, int *m, double *a,
 354         double *b, double *u, double *v, double *x, int *
 355         ix, int *job)
 356 {
 357     /* System generated locals */
 358     int i__1;
 359     double d__1;
 360
 361     /* Local variables */
 362     int i__, k;
 363
 364     /* Parameter adjustments */
 365     --ix;
 366     --x;
 367     --v;
 368     --u;
 369     --b;
 370     --a;
 371
 372     /* Function Body */
 373     k = 1;
 374     i__1 = *m;
 375     for (i__ = 1; i__ <= i__1; ++i__) {
 376         v[i__] = u[i__] * luksan_mxudot__(n, &x[1], &a[k], &ix[1], job);
 377         d__1 = -v[i__];
 378         luksan_mxudir__(n, &d__1, &b[k], &x[1], &x[1], &ix[1], job);
 379         k += *n;
 380 /* L1: */
 381     }
 382     return;
 383 } /* luksan_mxdrcb__ */
 384
 385 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 386 /* SUBROUTINE MXDRCF               ALL SYSTEMS                91/12/01
 387 * PURPOSE :
 388 * FORWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION OF
 389 * THE VECTOR X BY AN IMPLICIT BFGS UPDATE.
 390 *
 391 * PARAMETERS :
 392 *  II  N  NUMBER OF ROWS OF THE MATRICES A AND B.
 393 *  II  M  NUMBER OF COLUMNS OF THE MATRICES A AND B.
 394 *  RI  A(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 395 *  RI  B(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 396 *  RI  U(M)  VECTOR OF SCALAR COEFFICIENTS.
 397 *  RI  V(M)  VECTOR OF SCALAR COEFFICIENTS.
 398 *  RU  X(N)  PREMULTIPLIED VECTOR.
 399 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
 400 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
 401 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
 402 *         IX(I).EQ.-5.
 403 *
 404 * SUBPROGRAM USED :
 405 *  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
 406 *  RF  MXUDOT  DOT PRODUCT OF VECTORS.
 407 *
 408 * METHOD :
 409 * H.MATTHIES, G.STRANG: THE SOLUTION OF NONLINEAR FINITE ELEMENT
 410 * EQUATIONS. INT.J.NUMER. METHODS ENGN. 14 (1979) 1613-1626.
 411 */
 412 void luksan_mxdrcf__(int *n, int *m, double *a,
 413         double *b, double *u, double *v, double *x, int *
 414         ix, int *job)
 415 {
 416     /* System generated locals */
 417     double d__1;
 418
 419     /* Local variables */
 420     int i__, k;
 421     double temp;
 422
 423     /* Parameter adjustments */
 424     --ix;
 425     --x;
 426     --v;
 427     --u;
 428     --b;
 429     --a;
 430
 431     /* Function Body */
 432     k = (*m - 1) * *n + 1;
 433     for (i__ = *m; i__ >= 1; --i__) {
 434         temp = u[i__] * luksan_mxudot__(n, &x[1], &b[k], &ix[1], job);
 435         d__1 = v[i__] - temp;
 436         luksan_mxudir__(n, &d__1, &a[k], &x[1], &x[1], &ix[1], job);
 437         k -= *n;
 438 /* L1: */
 439     }
 440     return;
 441 } /* luksan_mxdrcf__ */
 442
 443 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 444 /* SUBROUTINE MXDRMM               ALL SYSTEMS                91/12/01
 445 * PURPOSE :
 446 * MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY
 447 * A VECTOR X.
 448 *
 449 * PARAMETERS :
 450 *  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
 451 *  II  M  NUMBER OF ROWS OF THE MATRIX A.
 452 *  RI  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
 453 *         ONE-DIMENSIONAL ARRAY.
 454 *  RI  X(N)  INPUT VECTOR.
 455 *  RO  Y(M)  OUTPUT VECTOR EQUAL TO A*X.
 456 */
 457 void luksan_mxdrmm__(int *n, int *m, double *a,
 458         double *x, double *y)
 459 {
 460     /* System generated locals */
 461     int i__1, i__2;
 462
 463     /* Local variables */
 464     int i__, j, k;
 465     double temp;
 466
 467     /* Parameter adjustments */
 468     --y;
 469     --x;
 470     --a;
 471
 472     /* Function Body */
 473     k = 0;
 474     i__1 = *m;
 475     for (j = 1; j <= i__1; ++j) {
 476         temp = 0.;
 477         i__2 = *n;
 478         for (i__ = 1; i__ <= i__2; ++i__) {
 479             temp += a[k + i__] * x[i__];
 480 /* L1: */
 481         }
 482         y[j] = temp;
 483         k += *n;
 484 /* L2: */
 485     }
 486     return;
 487 } /* luksan_mxdrmm__ */
 488
 489 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 490 /* SUBROUTINE MXDRSU               ALL SYSTEMS                91/12/01
 491 * PURPOSE :
 492 * SHIFT OF COLUMNS OF THE RECTANGULAR MATRICES A AND B. SHIFT OF
 493 * ELEMENTS OF THE VECTOR U. THESE SHIFTS ARE USED IN THE LIMITED
 494 * MEMORY BFGS METHOD.
 495 *
 496 * PARAMETERS :
 497 *  II  N  NUMBER OF ROWS OF THE MATRIX A AND B.
 498 *  II  M  NUMBER OF COLUMNS OF THE MATRIX A AND B.
 499 *  RU  A(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 500 *  RU  B(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 501 *  RU  U(M)  VECTOR.
 502 */
 503 void luksan_mxdrsu__(int *n, int *m, double *a,
 504                      double *b, double *u)
 505 {
 506     int i__, k, l;
 507
 508     /* Parameter adjustments */
 509     --u;
 510     --b;
 511     --a;
 512
 513     /* Function Body */
 514     k = (*m - 1) * *n + 1;
 515     for (i__ = *m - 1; i__ >= 1; --i__) {
 516         l = k - *n;
 517         luksan_mxvcop__(n, &a[l], &a[k]);
 518         luksan_mxvcop__(n, &b[l], &b[k]);
 519         u[i__ + 1] = u[i__];
 520         k = l;
 521 /* L1: */
 522     }
 523     return;
 524 } /* luksan_mxdrsu__ */
 525
 526 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 527 /* SUBROUTINE MXUCOP                ALL SYSTEMS                99/12/01
 528 * PURPOSE :
 529 * COPY OF THE VECTOR WITH INITIATION OF THE ACTIVE PART.
 530 *
 531 * PARAMETERS :
 532 *  II  N  VECTOR DIMENSION.
 533 *  RI  X(N)  INPUT VECTOR.
 534 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:= X.
 535 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
 536 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
 537 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
 538 *         IX(I).EQ.-5.
 539 */
 540 void luksan_mxucop__(int *n, double *x, double *y,
 541          int *ix, int *job)
 542 {
 543     /* System generated locals */
 544     int i__1;
 545
 546     /* Local variables */
 547     int i__;
 548
 549     /* Parameter adjustments */
 550     --ix;
 551     --y;
 552     --x;
 553
 554     /* Function Body */
 555     if (*job == 0) {
 556         i__1 = *n;
 557         for (i__ = 1; i__ <= i__1; ++i__) {
 558             y[i__] = x[i__];
 559 /* L1: */
 560         }
 561     } else if (*job > 0) {
 562         i__1 = *n;
 563         for (i__ = 1; i__ <= i__1; ++i__) {
 564             if (ix[i__] >= 0) {
 565                 y[i__] = x[i__];
 566             } else {
 567                 y[i__] = 0.;
 568             }
 569 /* L2: */
 570         }
 571     } else {
 572         i__1 = *n;
 573         for (i__ = 1; i__ <= i__1; ++i__) {
 574             if (ix[i__] != -5) {
 575                 y[i__] = x[i__];
 576             } else {
 577                 y[i__] = 0.;
 578             }
 579 /* L3: */
 580         }
 581     }
 582     return;
 583 } /* luksan_mxucop__ */
 584
 585 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 586 /* SUBROUTINE MXUDIR                ALL SYSTEMS                99/12/01
 587 * PURPOSE :
 588 * VECTOR AUGMENTED BY THE SCALED VECTOR IN A BOUND CONSTRAINED CASE.
 589 *
 590 * PARAMETERS :
 591 *  II  N  VECTOR DIMENSION.
 592 *  RI  A  SCALING FACTOR.
 593 *  RI  X(N)  INPUT VECTOR.
 594 *  RI  Y(N)  INPUT VECTOR.
 595 *  RO  Z(N)  OUTPUT VECTOR WHERE Z:= Y + A*X.
 596 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
 597 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
 598 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
 599 *         IX(I).EQ.-5.
 600 */
 601 void luksan_mxudir__(int *n, double *a, double *x,
 602          double *y, double *z__, int *ix, int *job)
 603 {
 604     /* System generated locals */
 605     int i__1;
 606
 607     /* Local variables */
 608     int i__;
 609
 610     /* Parameter adjustments */
 611     --ix;
 612     --z__;
 613     --y;
 614     --x;
 615
 616     /* Function Body */
 617     if (*job == 0) {
 618         i__1 = *n;
 619         for (i__ = 1; i__ <= i__1; ++i__) {
 620             z__[i__] = y[i__] + *a * x[i__];
 621 /* L1: */
 622         }
 623     } else if (*job > 0) {
 624         i__1 = *n;
 625         for (i__ = 1; i__ <= i__1; ++i__) {
 626             if (ix[i__] >= 0) {
 627                 z__[i__] = y[i__] + *a * x[i__];
 628             }
 629 /* L2: */
 630         }
 631     } else {
 632         i__1 = *n;
 633         for (i__ = 1; i__ <= i__1; ++i__) {
 634             if (ix[i__] != -5) {
 635                 z__[i__] = y[i__] + *a * x[i__];
 636             }
 637 /* L3: */
 638         }
 639     }
 640     return;
 641 } /* luksan_mxudir__ */
 642
 643 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 644 /* FUNCTION MXVDOT                  ALL SYSTEMS                91/12/01 */
 645 /* PURPOSE : */
 646 /* DOT PRODUCT OF TWO VECTORS. */
 647
 648 /* PARAMETERS : */
 649 /*  II  N  VECTOR DIMENSION. */
 650 /*  RI  X(N)  INPUT VECTOR. */
 651 /*  RI  Y(N)  INPUT VECTOR. */
 652 /*  RR  MXVDOT  VALUE OF DOT PRODUCT MXVDOT=TRANS(X)*Y. */
 653
 654 double luksan_mxvdot__(int *n, double *x, double *y)
 655 {
 656     /* System generated locals */
 657     int i__1;
 658
 659     /* Local variables */
 660     int i__;
 661     double temp;
 662
 663     /* Parameter adjustments */
 664     --y;
 665     --x;
 666
 667     /* Function Body */
 668     temp = 0.;
 669     i__1 = *n;
 670     for (i__ = 1; i__ <= i__1; ++i__) {
 671         temp += x[i__] * y[i__];
 672 /* L10: */
 673     }
 674     return temp;
 675 } /* luksan_mxvdot__ */
 676
 677 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 678 /* FUNCTION MXUDOT                  ALL SYSTEMS                99/12/01
 679 * PURPOSE :
 680 * DOT PRODUCT OF VECTORS IN A BOUND CONSTRAINED CASE.
 681 *
 682 * PARAMETERS :
 683 *  II  N  VECTOR DIMENSION.
 684 *  RI  X(N)  INPUT VECTOR.
 685 *  RI  Y(N)  INPUT VECTOR.
 686 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
 687 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
 688 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
 689 *         IX(I).EQ.-5.
 690 *  RR  MXUDOT  VALUE OF DOT PRODUCT MXUDOT=TRANS(X)*Y.
 691 */
 692 double luksan_mxudot__(int *n, double *x, double *y, int *ix,
 693                        int *job)
 694 {
 695     /* System generated locals */
 696     int i__1;
 697
 698     /* Local variables */
 699     int i__;
 700     double temp;
 701
 702     /* Parameter adjustments */
 703     --ix;
 704     --y;
 705     --x;
 706
 707     /* Function Body */
 708     temp = 0.;
 709     if (*job == 0) {
 710         i__1 = *n;
 711         for (i__ = 1; i__ <= i__1; ++i__) {
 712             temp += x[i__] * y[i__];
 713 /* L1: */
 714         }
 715     } else if (*job > 0) {
 716         i__1 = *n;
 717         for (i__ = 1; i__ <= i__1; ++i__) {
 718             if (ix[i__] >= 0) {
 719                 temp += x[i__] * y[i__];
 720             }
 721 /* L2: */
 722         }
 723     } else {
 724         i__1 = *n;
 725         for (i__ = 1; i__ <= i__1; ++i__) {
 726             if (ix[i__] != -5) {
 727                 temp += x[i__] * y[i__];
 728             }
 729 /* L3: */
 730         }
 731     }
 732     return temp;
 733 } /* luksan_mxudot__ */
 734
 735 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 736 /* SUBROUTINE MXUNEG                ALL SYSTEMS                00/12/01
 737 * PURPOSE :
 738 * COPY OF THE VECTOR WITH INITIATION OF THE ACTIVE PART.
 739 *
 740 * PARAMETERS :
 741 *  II  N  VECTOR DIMENSION.
 742 *  RI  X(N)  INPUT VECTOR.
 743 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:= X.
 744 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
 745 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
 746 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
 747 *         IX(I).EQ.-5.
 748 */
 749 void luksan_mxuneg__(int *n, double *x, double *y,
 750          int *ix, int *job)
 751 {
 752     /* System generated locals */
 753     int i__1;
 754
 755     /* Local variables */
 756     int i__;
 757
 758     /* Parameter adjustments */
 759     --ix;
 760     --y;
 761     --x;
 762
 763     /* Function Body */
 764     if (*job == 0) {
 765         i__1 = *n;
 766         for (i__ = 1; i__ <= i__1; ++i__) {
 767             y[i__] = -x[i__];
 768 /* L1: */
 769         }
 770     } else if (*job > 0) {
 771         i__1 = *n;
 772         for (i__ = 1; i__ <= i__1; ++i__) {
 773             if (ix[i__] >= 0) {
 774                 y[i__] = -x[i__];
 775             } else {
 776                 y[i__] = 0.;
 777             }
 778 /* L2: */
 779         }
 780     } else {
 781         i__1 = *n;
 782         for (i__ = 1; i__ <= i__1; ++i__) {
 783             if (ix[i__] != -5) {
 784                 y[i__] = -x[i__];
 785             } else {
 786                 y[i__] = 0.;
 787             }
 788 /* L3: */
 789         }
 790     }
 791     return;
 792 } /* luksan_mxuneg__ */
 793
 794 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 795 /* SUBROUTINE MXUZER                ALL SYSTEMS                99/12/01
 796 * PURPOSE :
 797 * VECTOR ELEMENTS CORRESPONDING TO ACTIVE BOUNDS ARE SET TO ZERO.
 798 *
 799 * PARAMETERS :
 800 *  II  N  VECTOR DIMENSION.
 801 *  RO  X(N)  OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I.
 802 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
 803 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
 804 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
 805 *         IX(I).EQ.-5.
 806 */
 807 void luksan_mxuzer__(int *n, double *x, int *ix,
 808                      int *job)
 809 {
 810     /* System generated locals */
 811     int i__1;
 812
 813     /* Local variables */
 814     int i__;
 815
 816     /* Parameter adjustments */
 817     --ix;
 818     --x;
 819
 820     /* Function Body */
 821     if (*job == 0) {
 822         return;
 823     }
 824     i__1 = *n;
 825     for (i__ = 1; i__ <= i__1; ++i__) {
 826         if (ix[i__] < 0) {
 827             x[i__] = 0.;
 828         }
 829 /* L1: */
 830     }
 831     return;
 832 } /* luksan_mxuzer__ */
 833
 834 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 835 /* SUBROUTINE MXVCOP                ALL SYSTEMS                88/12/01
 836 * PURPOSE :
 837 * COPYING OF A VECTOR.
 838 *
 839 * PARAMETERS :
 840 *  II  N  VECTOR DIMENSION.
 841 *  RI  X(N)  INPUT VECTOR.
 842 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:= X.
 843 */
 844 void luksan_mxvcop__(int *n, double *x, double *y)
 845 {
 846     /* System generated locals */
 847     int i__1;
 848
 849     /* Local variables */
 850     int i__;
 851
 852     /* Parameter adjustments */
 853     --y;
 854     --x;
 855
 856     /* Function Body */
 857     i__1 = *n;
 858     for (i__ = 1; i__ <= i__1; ++i__) {
 859         y[i__] = x[i__];
 860 /* L10: */
 861     }
 862     return;
 863 } /* luksan_mxvcop__ */
 864
 865 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 866 /* SUBROUTINE MXVDIF                ALL SYSTEMS                88/12/01
 867 * PURPOSE :
 868 * VECTOR DIFFERENCE.
 869 *
 870 * PARAMETERS :
 871 *  RI  X(N)  INPUT VECTOR.
 872 *  RI  Y(N)  INPUT VECTOR.
 873 *  RO  Z(N)  OUTPUT VECTOR WHERE Z:= X - Y.
 874 */
 875 void luksan_mxvdif__(int *n, double *x, double *y,
 876          double *z__)
 877 {
 878     /* System generated locals */
 879     int i__1;
 880
 881     /* Local variables */
 882     int i__;
 883
 884     /* Parameter adjustments */
 885     --z__;
 886     --y;
 887     --x;
 888
 889     /* Function Body */
 890     i__1 = *n;
 891     for (i__ = 1; i__ <= i__1; ++i__) {
 892         z__[i__] = x[i__] - y[i__];
 893 /* L10: */
 894     }
 895     return;
 896 } /* luksan_mxvdif__ */
 897
 898 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 899 /* SUBROUTINE MXVNEG                ALL SYSTEMS                88/12/01
 900 * PURPOSE :
 901 * CHANGE THE SIGNS OF VECTOR ELEMENTS.
 902 *
 903 * PARAMETERS :
 904 *  II  N  VECTOR DIMENSION.
 905 *  RI  X(N)  INPUT VECTOR.
 906 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:= - X.
 907 */
 908 void luksan_mxvneg__(int *n, double *x, double *y)
 909 {
 910     /* System generated locals */
 911     int i__1;
 912
 913     /* Local variables */
 914     int i__;
 915
 916     /* Parameter adjustments */
 917     --y;
 918     --x;
 919
 920     /* Function Body */
 921     i__1 = *n;
 922     for (i__ = 1; i__ <= i__1; ++i__) {
 923         y[i__] = -x[i__];
 924 /* L10: */
 925     }
 926     return;
 927 } /* luksan_mxvneg__ */
 928
 929 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 930 /* SUBROUTINE MXVSCL                ALL SYSTEMS                88/12/01
 931 * PURPOSE :
 932 * SCALING OF A VECTOR.
 933 *
 934 * PARAMETERS :
 935 *  II  N  VECTOR DIMENSION.
 936 *  RI  X(N)  INPUT VECTOR.
 937 *  RI  A  SCALING FACTOR.
 938 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:= A*X.
 939 */
 940 void luksan_mxvscl__(int *n, double *a, double *x,
 941          double *y)
 942 {
 943     /* System generated locals */
 944     int i__1;
 945
 946     /* Local variables */
 947     int i__;
 948
 949     /* Parameter adjustments */
 950     --y;
 951     --x;
 952
 953     /* Function Body */
 954     i__1 = *n;
 955     for (i__ = 1; i__ <= i__1; ++i__) {
 956         y[i__] = *a * x[i__];
 957 /* L1: */
 958     }
 959     return;
 960 } /* luksan_mxvscl__ */
 961
 962 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
 963 /* SUBROUTINE MXVSET                ALL SYSTEMS                88/12/01
 964 * PURPOSE :
 965 * A SCALAR IS SET TO ALL THE ELEMENTS OF A VECTOR.
 966 *
 967 * PARAMETERS :
 968 *  II  N  VECTOR DIMENSION.
 969 *  RI  A  INITIAL VALUE.
 970 *  RO  X(N)  OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I.
 971 */
 972 void luksan_mxvset__(int *n, double *a, double *x)
 973 {
 974     /* System generated locals */
 975     int i__1;
 976
 977     /* Local variables */
 978     int i__;
 979
 980     /* Parameter adjustments */
 981     --x;
 982
 983     /* Function Body */
 984     i__1 = *n;
 985     for (i__ = 1; i__ <= i__1; ++i__) {
 986         x[i__] = *a;
 987 /* L10: */
 988     }
 989     return;
 990 } /* luksan_mxvset__ */
 991