slansf.c
Go to the documentation of this file.
00001 /* slansf.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__1 = 1;
00019 
00020 doublereal slansf_(char *norm, char *transr, char *uplo, integer *n, real *a, 
00021         real *work)
00022 {
00023     /* System generated locals */
00024     integer i__1, i__2;
00025     real ret_val, r__1, r__2, r__3;
00026 
00027     /* Builtin functions */
00028     double sqrt(doublereal);
00029 
00030     /* Local variables */
00031     integer i__, j, k, l;
00032     real s;
00033     integer n1;
00034     real aa;
00035     integer lda, ifm, noe, ilu;
00036     real scale;
00037     extern logical lsame_(char *, char *);
00038     real value;
00039     extern integer isamax_(integer *, real *, integer *);
00040     extern /* Subroutine */ int slassq_(integer *, real *, integer *, real *, 
00041             real *);
00042 
00043 
00044 /*  -- LAPACK routine (version 3.2)                                    -- */
00045 
00046 /*  -- Contributed by Fred Gustavson of the IBM Watson Research Center -- */
00047 /*  -- November 2008                                                   -- */
00048 
00049 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00050 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00051 
00052 /*     .. Scalar Arguments .. */
00053 /*     .. */
00054 /*     .. Array Arguments .. */
00055 /*     .. */
00056 
00057 /*  Purpose */
00058 /*  ======= */
00059 
00060 /*  SLANSF returns the value of the one norm, or the Frobenius norm, or */
00061 /*  the infinity norm, or the element of largest absolute value of a */
00062 /*  real symmetric matrix A in RFP format. */
00063 
00064 /*  Description */
00065 /*  =========== */
00066 
00067 /*  SLANSF returns the value */
00068 
00069 /*     SLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm' */
00070 /*              ( */
00071 /*              ( norm1(A),         NORM = '1', 'O' or 'o' */
00072 /*              ( */
00073 /*              ( normI(A),         NORM = 'I' or 'i' */
00074 /*              ( */
00075 /*              ( normF(A),         NORM = 'F', 'f', 'E' or 'e' */
00076 
00077 /*  where  norm1  denotes the  one norm of a matrix (maximum column sum), */
00078 /*  normI  denotes the  infinity norm  of a matrix  (maximum row sum) and */
00079 /*  normF  denotes the  Frobenius norm of a matrix (square root of sum of */
00080 /*  squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm. */
00081 
00082 /*  Arguments */
00083 /*  ========= */
00084 
00085 /*  NORM    (input) CHARACTER */
00086 /*          Specifies the value to be returned in SLANSF as described */
00087 /*          above. */
00088 
00089 /*  TRANSR  (input) CHARACTER */
00090 /*          Specifies whether the RFP format of A is normal or */
00091 /*          transposed format. */
00092 /*          = 'N':  RFP format is Normal; */
00093 /*          = 'T':  RFP format is Transpose. */
00094 
00095 /*  UPLO    (input) CHARACTER */
00096 /*           On entry, UPLO specifies whether the RFP matrix A came from */
00097 /*           an upper or lower triangular matrix as follows: */
00098 /*           = 'U': RFP A came from an upper triangular matrix; */
00099 /*           = 'L': RFP A came from a lower triangular matrix. */
00100 
00101 /*  N       (input) INTEGER */
00102 /*          The order of the matrix A. N >= 0. When N = 0, SLANSF is */
00103 /*          set to zero. */
00104 
00105 /*  A       (input) REAL array, dimension ( N*(N+1)/2 ); */
00106 /*          On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L') */
00107 /*          part of the symmetric matrix A stored in RFP format. See the */
00108 /*          "Notes" below for more details. */
00109 /*          Unchanged on exit. */
00110 
00111 /*  WORK    (workspace) REAL array, dimension (MAX(1,LWORK)), */
00112 /*          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, */
00113 /*          WORK is not referenced. */
00114 
00115 /*  Notes */
00116 /*  ===== */
00117 
00118 /*  We first consider Rectangular Full Packed (RFP) Format when N is */
00119 /*  even. We give an example where N = 6. */
00120 
00121 /*      AP is Upper             AP is Lower */
00122 
00123 /*   00 01 02 03 04 05       00 */
00124 /*      11 12 13 14 15       10 11 */
00125 /*         22 23 24 25       20 21 22 */
00126 /*            33 34 35       30 31 32 33 */
00127 /*               44 45       40 41 42 43 44 */
00128 /*                  55       50 51 52 53 54 55 */
00129 
00130 
00131 /*  Let TRANSR = 'N'. RFP holds AP as follows: */
00132 /*  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last */
00133 /*  three columns of AP upper. The lower triangle A(4:6,0:2) consists of */
00134 /*  the transpose of the first three columns of AP upper. */
00135 /*  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first */
00136 /*  three columns of AP lower. The upper triangle A(0:2,0:2) consists of */
00137 /*  the transpose of the last three columns of AP lower. */
00138 /*  This covers the case N even and TRANSR = 'N'. */
00139 
00140 /*         RFP A                   RFP A */
00141 
00142 /*        03 04 05                33 43 53 */
00143 /*        13 14 15                00 44 54 */
00144 /*        23 24 25                10 11 55 */
00145 /*        33 34 35                20 21 22 */
00146 /*        00 44 45                30 31 32 */
00147 /*        01 11 55                40 41 42 */
00148 /*        02 12 22                50 51 52 */
00149 
00150 /*  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the */
00151 /*  transpose of RFP A above. One therefore gets: */
00152 
00153 
00154 /*           RFP A                   RFP A */
00155 
00156 /*     03 13 23 33 00 01 02    33 00 10 20 30 40 50 */
00157 /*     04 14 24 34 44 11 12    43 44 11 21 31 41 51 */
00158 /*     05 15 25 35 45 55 22    53 54 55 22 32 42 52 */
00159 
00160 
00161 /*  We first consider Rectangular Full Packed (RFP) Format when N is */
00162 /*  odd. We give an example where N = 5. */
00163 
00164 /*     AP is Upper                 AP is Lower */
00165 
00166 /*   00 01 02 03 04              00 */
00167 /*      11 12 13 14              10 11 */
00168 /*         22 23 24              20 21 22 */
00169 /*            33 34              30 31 32 33 */
00170 /*               44              40 41 42 43 44 */
00171 
00172 
00173 /*  Let TRANSR = 'N'. RFP holds AP as follows: */
00174 /*  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last */
00175 /*  three columns of AP upper. The lower triangle A(3:4,0:1) consists of */
00176 /*  the transpose of the first two columns of AP upper. */
00177 /*  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first */
00178 /*  three columns of AP lower. The upper triangle A(0:1,1:2) consists of */
00179 /*  the transpose of the last two columns of AP lower. */
00180 /*  This covers the case N odd and TRANSR = 'N'. */
00181 
00182 /*         RFP A                   RFP A */
00183 
00184 /*        02 03 04                00 33 43 */
00185 /*        12 13 14                10 11 44 */
00186 /*        22 23 24                20 21 22 */
00187 /*        00 33 34                30 31 32 */
00188 /*        01 11 44                40 41 42 */
00189 
00190 /*  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the */
00191 /*  transpose of RFP A above. One therefore gets: */
00192 
00193 /*           RFP A                   RFP A */
00194 
00195 /*     02 12 22 00 01             00 10 20 30 40 50 */
00196 /*     03 13 23 33 11             33 11 21 31 41 51 */
00197 /*     04 14 24 34 44             43 44 22 32 42 52 */
00198 
00199 /*  Reference */
00200 /*  ========= */
00201 
00202 /*  ===================================================================== */
00203 
00204 /*     .. */
00205 /*     .. Parameters .. */
00206 /*     .. */
00207 /*     .. Local Scalars .. */
00208 /*     .. */
00209 /*     .. External Functions .. */
00210 /*     .. */
00211 /*     .. External Subroutines .. */
00212 /*     .. */
00213 /*     .. Intrinsic Functions .. */
00214 /*     .. */
00215 /*     .. Executable Statements .. */
00216 
00217     if (*n == 0) {
00218         ret_val = 0.f;
00219         return ret_val;
00220     }
00221 
00222 /*     set noe = 1 if n is odd. if n is even set noe=0 */
00223 
00224     noe = 1;
00225     if (*n % 2 == 0) {
00226         noe = 0;
00227     }
00228 
00229 /*     set ifm = 0 when form='T or 't' and 1 otherwise */
00230 
00231     ifm = 1;
00232     if (lsame_(transr, "T")) {
00233         ifm = 0;
00234     }
00235 
00236 /*     set ilu = 0 when uplo='U or 'u' and 1 otherwise */
00237 
00238     ilu = 1;
00239     if (lsame_(uplo, "U")) {
00240         ilu = 0;
00241     }
00242 
00243 /*     set lda = (n+1)/2 when ifm = 0 */
00244 /*     set lda = n when ifm = 1 and noe = 1 */
00245 /*     set lda = n+1 when ifm = 1 and noe = 0 */
00246 
00247     if (ifm == 1) {
00248         if (noe == 1) {
00249             lda = *n;
00250         } else {
00251 /*           noe=0 */
00252             lda = *n + 1;
00253         }
00254     } else {
00255 /*        ifm=0 */
00256         lda = (*n + 1) / 2;
00257     }
00258 
00259     if (lsame_(norm, "M")) {
00260 
00261 /*       Find max(abs(A(i,j))). */
00262 
00263         k = (*n + 1) / 2;
00264         value = 0.f;
00265         if (noe == 1) {
00266 /*           n is odd */
00267             if (ifm == 1) {
00268 /*           A is n by k */
00269                 i__1 = k - 1;
00270                 for (j = 0; j <= i__1; ++j) {
00271                     i__2 = *n - 1;
00272                     for (i__ = 0; i__ <= i__2; ++i__) {
00273 /* Computing MAX */
00274                         r__2 = value, r__3 = (r__1 = a[i__ + j * lda], dabs(
00275                                 r__1));
00276                         value = dmax(r__2,r__3);
00277                     }
00278                 }
00279             } else {
00280 /*              xpose case; A is k by n */
00281                 i__1 = *n - 1;
00282                 for (j = 0; j <= i__1; ++j) {
00283                     i__2 = k - 1;
00284                     for (i__ = 0; i__ <= i__2; ++i__) {
00285 /* Computing MAX */
00286                         r__2 = value, r__3 = (r__1 = a[i__ + j * lda], dabs(
00287                                 r__1));
00288                         value = dmax(r__2,r__3);
00289                     }
00290                 }
00291             }
00292         } else {
00293 /*           n is even */
00294             if (ifm == 1) {
00295 /*              A is n+1 by k */
00296                 i__1 = k - 1;
00297                 for (j = 0; j <= i__1; ++j) {
00298                     i__2 = *n;
00299                     for (i__ = 0; i__ <= i__2; ++i__) {
00300 /* Computing MAX */
00301                         r__2 = value, r__3 = (r__1 = a[i__ + j * lda], dabs(
00302                                 r__1));
00303                         value = dmax(r__2,r__3);
00304                     }
00305                 }
00306             } else {
00307 /*              xpose case; A is k by n+1 */
00308                 i__1 = *n;
00309                 for (j = 0; j <= i__1; ++j) {
00310                     i__2 = k - 1;
00311                     for (i__ = 0; i__ <= i__2; ++i__) {
00312 /* Computing MAX */
00313                         r__2 = value, r__3 = (r__1 = a[i__ + j * lda], dabs(
00314                                 r__1));
00315                         value = dmax(r__2,r__3);
00316                     }
00317                 }
00318             }
00319         }
00320     } else if (lsame_(norm, "I") || lsame_(norm, "O") || *(unsigned char *)norm == '1') {
00321 
00322 /*        Find normI(A) ( = norm1(A), since A is symmetric). */
00323 
00324         if (ifm == 1) {
00325             k = *n / 2;
00326             if (noe == 1) {
00327 /*              n is odd */
00328                 if (ilu == 0) {
00329                     i__1 = k - 1;
00330                     for (i__ = 0; i__ <= i__1; ++i__) {
00331                         work[i__] = 0.f;
00332                     }
00333                     i__1 = k;
00334                     for (j = 0; j <= i__1; ++j) {
00335                         s = 0.f;
00336                         i__2 = k + j - 1;
00337                         for (i__ = 0; i__ <= i__2; ++i__) {
00338                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00339 /*                       -> A(i,j+k) */
00340                             s += aa;
00341                             work[i__] += aa;
00342                         }
00343                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00344 /*                    -> A(j+k,j+k) */
00345                         work[j + k] = s + aa;
00346                         if (i__ == k + k) {
00347                             goto L10;
00348                         }
00349                         ++i__;
00350                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00351 /*                    -> A(j,j) */
00352                         work[j] += aa;
00353                         s = 0.f;
00354                         i__2 = k - 1;
00355                         for (l = j + 1; l <= i__2; ++l) {
00356                             ++i__;
00357                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00358 /*                       -> A(l,j) */
00359                             s += aa;
00360                             work[l] += aa;
00361                         }
00362                         work[j] += s;
00363                     }
00364 L10:
00365                     i__ = isamax_(n, work, &c__1);
00366                     value = work[i__ - 1];
00367                 } else {
00368 /*                 ilu = 1 */
00369                     ++k;
00370 /*                 k=(n+1)/2 for n odd and ilu=1 */
00371                     i__1 = *n - 1;
00372                     for (i__ = k; i__ <= i__1; ++i__) {
00373                         work[i__] = 0.f;
00374                     }
00375                     for (j = k - 1; j >= 0; --j) {
00376                         s = 0.f;
00377                         i__1 = j - 2;
00378                         for (i__ = 0; i__ <= i__1; ++i__) {
00379                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00380 /*                       -> A(j+k,i+k) */
00381                             s += aa;
00382                             work[i__ + k] += aa;
00383                         }
00384                         if (j > 0) {
00385                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00386 /*                       -> A(j+k,j+k) */
00387                             s += aa;
00388                             work[i__ + k] += s;
00389 /*                       i=j */
00390                             ++i__;
00391                         }
00392                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00393 /*                    -> A(j,j) */
00394                         work[j] = aa;
00395                         s = 0.f;
00396                         i__1 = *n - 1;
00397                         for (l = j + 1; l <= i__1; ++l) {
00398                             ++i__;
00399                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00400 /*                       -> A(l,j) */
00401                             s += aa;
00402                             work[l] += aa;
00403                         }
00404                         work[j] += s;
00405                     }
00406                     i__ = isamax_(n, work, &c__1);
00407                     value = work[i__ - 1];
00408                 }
00409             } else {
00410 /*              n is even */
00411                 if (ilu == 0) {
00412                     i__1 = k - 1;
00413                     for (i__ = 0; i__ <= i__1; ++i__) {
00414                         work[i__] = 0.f;
00415                     }
00416                     i__1 = k - 1;
00417                     for (j = 0; j <= i__1; ++j) {
00418                         s = 0.f;
00419                         i__2 = k + j - 1;
00420                         for (i__ = 0; i__ <= i__2; ++i__) {
00421                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00422 /*                       -> A(i,j+k) */
00423                             s += aa;
00424                             work[i__] += aa;
00425                         }
00426                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00427 /*                    -> A(j+k,j+k) */
00428                         work[j + k] = s + aa;
00429                         ++i__;
00430                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00431 /*                    -> A(j,j) */
00432                         work[j] += aa;
00433                         s = 0.f;
00434                         i__2 = k - 1;
00435                         for (l = j + 1; l <= i__2; ++l) {
00436                             ++i__;
00437                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00438 /*                       -> A(l,j) */
00439                             s += aa;
00440                             work[l] += aa;
00441                         }
00442                         work[j] += s;
00443                     }
00444                     i__ = isamax_(n, work, &c__1);
00445                     value = work[i__ - 1];
00446                 } else {
00447 /*                 ilu = 1 */
00448                     i__1 = *n - 1;
00449                     for (i__ = k; i__ <= i__1; ++i__) {
00450                         work[i__] = 0.f;
00451                     }
00452                     for (j = k - 1; j >= 0; --j) {
00453                         s = 0.f;
00454                         i__1 = j - 1;
00455                         for (i__ = 0; i__ <= i__1; ++i__) {
00456                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00457 /*                       -> A(j+k,i+k) */
00458                             s += aa;
00459                             work[i__ + k] += aa;
00460                         }
00461                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00462 /*                    -> A(j+k,j+k) */
00463                         s += aa;
00464                         work[i__ + k] += s;
00465 /*                    i=j */
00466                         ++i__;
00467                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00468 /*                    -> A(j,j) */
00469                         work[j] = aa;
00470                         s = 0.f;
00471                         i__1 = *n - 1;
00472                         for (l = j + 1; l <= i__1; ++l) {
00473                             ++i__;
00474                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00475 /*                       -> A(l,j) */
00476                             s += aa;
00477                             work[l] += aa;
00478                         }
00479                         work[j] += s;
00480                     }
00481                     i__ = isamax_(n, work, &c__1);
00482                     value = work[i__ - 1];
00483                 }
00484             }
00485         } else {
00486 /*           ifm=0 */
00487             k = *n / 2;
00488             if (noe == 1) {
00489 /*              n is odd */
00490                 if (ilu == 0) {
00491                     n1 = k;
00492 /*                 n/2 */
00493                     ++k;
00494 /*                 k is the row size and lda */
00495                     i__1 = *n - 1;
00496                     for (i__ = n1; i__ <= i__1; ++i__) {
00497                         work[i__] = 0.f;
00498                     }
00499                     i__1 = n1 - 1;
00500                     for (j = 0; j <= i__1; ++j) {
00501                         s = 0.f;
00502                         i__2 = k - 1;
00503                         for (i__ = 0; i__ <= i__2; ++i__) {
00504                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00505 /*                       A(j,n1+i) */
00506                             work[i__ + n1] += aa;
00507                             s += aa;
00508                         }
00509                         work[j] = s;
00510                     }
00511 /*                 j=n1=k-1 is special */
00512                     s = (r__1 = a[j * lda], dabs(r__1));
00513 /*                 A(k-1,k-1) */
00514                     i__1 = k - 1;
00515                     for (i__ = 1; i__ <= i__1; ++i__) {
00516                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00517 /*                    A(k-1,i+n1) */
00518                         work[i__ + n1] += aa;
00519                         s += aa;
00520                     }
00521                     work[j] += s;
00522                     i__1 = *n - 1;
00523                     for (j = k; j <= i__1; ++j) {
00524                         s = 0.f;
00525                         i__2 = j - k - 1;
00526                         for (i__ = 0; i__ <= i__2; ++i__) {
00527                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00528 /*                       A(i,j-k) */
00529                             work[i__] += aa;
00530                             s += aa;
00531                         }
00532 /*                    i=j-k */
00533                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00534 /*                    A(j-k,j-k) */
00535                         s += aa;
00536                         work[j - k] += s;
00537                         ++i__;
00538                         s = (r__1 = a[i__ + j * lda], dabs(r__1));
00539 /*                    A(j,j) */
00540                         i__2 = *n - 1;
00541                         for (l = j + 1; l <= i__2; ++l) {
00542                             ++i__;
00543                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00544 /*                       A(j,l) */
00545                             work[l] += aa;
00546                             s += aa;
00547                         }
00548                         work[j] += s;
00549                     }
00550                     i__ = isamax_(n, work, &c__1);
00551                     value = work[i__ - 1];
00552                 } else {
00553 /*                 ilu=1 */
00554                     ++k;
00555 /*                 k=(n+1)/2 for n odd and ilu=1 */
00556                     i__1 = *n - 1;
00557                     for (i__ = k; i__ <= i__1; ++i__) {
00558                         work[i__] = 0.f;
00559                     }
00560                     i__1 = k - 2;
00561                     for (j = 0; j <= i__1; ++j) {
00562 /*                    process */
00563                         s = 0.f;
00564                         i__2 = j - 1;
00565                         for (i__ = 0; i__ <= i__2; ++i__) {
00566                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00567 /*                       A(j,i) */
00568                             work[i__] += aa;
00569                             s += aa;
00570                         }
00571                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00572 /*                    i=j so process of A(j,j) */
00573                         s += aa;
00574                         work[j] = s;
00575 /*                    is initialised here */
00576                         ++i__;
00577 /*                    i=j process A(j+k,j+k) */
00578                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00579                         s = aa;
00580                         i__2 = *n - 1;
00581                         for (l = k + j + 1; l <= i__2; ++l) {
00582                             ++i__;
00583                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00584 /*                       A(l,k+j) */
00585                             s += aa;
00586                             work[l] += aa;
00587                         }
00588                         work[k + j] += s;
00589                     }
00590 /*                 j=k-1 is special :process col A(k-1,0:k-1) */
00591                     s = 0.f;
00592                     i__1 = k - 2;
00593                     for (i__ = 0; i__ <= i__1; ++i__) {
00594                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00595 /*                    A(k,i) */
00596                         work[i__] += aa;
00597                         s += aa;
00598                     }
00599 /*                 i=k-1 */
00600                     aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00601 /*                 A(k-1,k-1) */
00602                     s += aa;
00603                     work[i__] = s;
00604 /*                 done with col j=k+1 */
00605                     i__1 = *n - 1;
00606                     for (j = k; j <= i__1; ++j) {
00607 /*                    process col j of A = A(j,0:k-1) */
00608                         s = 0.f;
00609                         i__2 = k - 1;
00610                         for (i__ = 0; i__ <= i__2; ++i__) {
00611                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00612 /*                       A(j,i) */
00613                             work[i__] += aa;
00614                             s += aa;
00615                         }
00616                         work[j] += s;
00617                     }
00618                     i__ = isamax_(n, work, &c__1);
00619                     value = work[i__ - 1];
00620                 }
00621             } else {
00622 /*              n is even */
00623                 if (ilu == 0) {
00624                     i__1 = *n - 1;
00625                     for (i__ = k; i__ <= i__1; ++i__) {
00626                         work[i__] = 0.f;
00627                     }
00628                     i__1 = k - 1;
00629                     for (j = 0; j <= i__1; ++j) {
00630                         s = 0.f;
00631                         i__2 = k - 1;
00632                         for (i__ = 0; i__ <= i__2; ++i__) {
00633                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00634 /*                       A(j,i+k) */
00635                             work[i__ + k] += aa;
00636                             s += aa;
00637                         }
00638                         work[j] = s;
00639                     }
00640 /*                 j=k */
00641                     aa = (r__1 = a[j * lda], dabs(r__1));
00642 /*                 A(k,k) */
00643                     s = aa;
00644                     i__1 = k - 1;
00645                     for (i__ = 1; i__ <= i__1; ++i__) {
00646                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00647 /*                    A(k,k+i) */
00648                         work[i__ + k] += aa;
00649                         s += aa;
00650                     }
00651                     work[j] += s;
00652                     i__1 = *n - 1;
00653                     for (j = k + 1; j <= i__1; ++j) {
00654                         s = 0.f;
00655                         i__2 = j - 2 - k;
00656                         for (i__ = 0; i__ <= i__2; ++i__) {
00657                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00658 /*                       A(i,j-k-1) */
00659                             work[i__] += aa;
00660                             s += aa;
00661                         }
00662 /*                     i=j-1-k */
00663                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00664 /*                    A(j-k-1,j-k-1) */
00665                         s += aa;
00666                         work[j - k - 1] += s;
00667                         ++i__;
00668                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00669 /*                    A(j,j) */
00670                         s = aa;
00671                         i__2 = *n - 1;
00672                         for (l = j + 1; l <= i__2; ++l) {
00673                             ++i__;
00674                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00675 /*                       A(j,l) */
00676                             work[l] += aa;
00677                             s += aa;
00678                         }
00679                         work[j] += s;
00680                     }
00681 /*                 j=n */
00682                     s = 0.f;
00683                     i__1 = k - 2;
00684                     for (i__ = 0; i__ <= i__1; ++i__) {
00685                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00686 /*                    A(i,k-1) */
00687                         work[i__] += aa;
00688                         s += aa;
00689                     }
00690 /*                 i=k-1 */
00691                     aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00692 /*                 A(k-1,k-1) */
00693                     s += aa;
00694                     work[i__] += s;
00695                     i__ = isamax_(n, work, &c__1);
00696                     value = work[i__ - 1];
00697                 } else {
00698 /*                 ilu=1 */
00699                     i__1 = *n - 1;
00700                     for (i__ = k; i__ <= i__1; ++i__) {
00701                         work[i__] = 0.f;
00702                     }
00703 /*                 j=0 is special :process col A(k:n-1,k) */
00704                     s = dabs(a[0]);
00705 /*                 A(k,k) */
00706                     i__1 = k - 1;
00707                     for (i__ = 1; i__ <= i__1; ++i__) {
00708                         aa = (r__1 = a[i__], dabs(r__1));
00709 /*                    A(k+i,k) */
00710                         work[i__ + k] += aa;
00711                         s += aa;
00712                     }
00713                     work[k] += s;
00714                     i__1 = k - 1;
00715                     for (j = 1; j <= i__1; ++j) {
00716 /*                    process */
00717                         s = 0.f;
00718                         i__2 = j - 2;
00719                         for (i__ = 0; i__ <= i__2; ++i__) {
00720                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00721 /*                       A(j-1,i) */
00722                             work[i__] += aa;
00723                             s += aa;
00724                         }
00725                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00726 /*                    i=j-1 so process of A(j-1,j-1) */
00727                         s += aa;
00728                         work[j - 1] = s;
00729 /*                    is initialised here */
00730                         ++i__;
00731 /*                    i=j process A(j+k,j+k) */
00732                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00733                         s = aa;
00734                         i__2 = *n - 1;
00735                         for (l = k + j + 1; l <= i__2; ++l) {
00736                             ++i__;
00737                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00738 /*                       A(l,k+j) */
00739                             s += aa;
00740                             work[l] += aa;
00741                         }
00742                         work[k + j] += s;
00743                     }
00744 /*                 j=k is special :process col A(k,0:k-1) */
00745                     s = 0.f;
00746                     i__1 = k - 2;
00747                     for (i__ = 0; i__ <= i__1; ++i__) {
00748                         aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00749 /*                    A(k,i) */
00750                         work[i__] += aa;
00751                         s += aa;
00752                     }
00753 /*                 i=k-1 */
00754                     aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00755 /*                 A(k-1,k-1) */
00756                     s += aa;
00757                     work[i__] = s;
00758 /*                 done with col j=k+1 */
00759                     i__1 = *n;
00760                     for (j = k + 1; j <= i__1; ++j) {
00761 /*                    process col j-1 of A = A(j-1,0:k-1) */
00762                         s = 0.f;
00763                         i__2 = k - 1;
00764                         for (i__ = 0; i__ <= i__2; ++i__) {
00765                             aa = (r__1 = a[i__ + j * lda], dabs(r__1));
00766 /*                       A(j-1,i) */
00767                             work[i__] += aa;
00768                             s += aa;
00769                         }
00770                         work[j - 1] += s;
00771                     }
00772                     i__ = isamax_(n, work, &c__1);
00773                     value = work[i__ - 1];
00774                 }
00775             }
00776         }
00777     } else if (lsame_(norm, "F") || lsame_(norm, "E")) {
00778 
00779 /*       Find normF(A). */
00780 
00781         k = (*n + 1) / 2;
00782         scale = 0.f;
00783         s = 1.f;
00784         if (noe == 1) {
00785 /*           n is odd */
00786             if (ifm == 1) {
00787 /*              A is normal */
00788                 if (ilu == 0) {
00789 /*                 A is upper */
00790                     i__1 = k - 3;
00791                     for (j = 0; j <= i__1; ++j) {
00792                         i__2 = k - j - 2;
00793                         slassq_(&i__2, &a[k + j + 1 + j * lda], &c__1, &scale, 
00794                                  &s);
00795 /*                    L at A(k,0) */
00796                     }
00797                     i__1 = k - 1;
00798                     for (j = 0; j <= i__1; ++j) {
00799                         i__2 = k + j - 1;
00800                         slassq_(&i__2, &a[j * lda], &c__1, &scale, &s);
00801 /*                    trap U at A(0,0) */
00802                     }
00803                     s += s;
00804 /*                 double s for the off diagonal elements */
00805                     i__1 = k - 1;
00806                     i__2 = lda + 1;
00807                     slassq_(&i__1, &a[k], &i__2, &scale, &s);
00808 /*                 tri L at A(k,0) */
00809                     i__1 = lda + 1;
00810                     slassq_(&k, &a[k - 1], &i__1, &scale, &s);
00811 /*                 tri U at A(k-1,0) */
00812                 } else {
00813 /*                 ilu=1 & A is lower */
00814                     i__1 = k - 1;
00815                     for (j = 0; j <= i__1; ++j) {
00816                         i__2 = *n - j - 1;
00817                         slassq_(&i__2, &a[j + 1 + j * lda], &c__1, &scale, &s)
00818                                 ;
00819 /*                    trap L at A(0,0) */
00820                     }
00821                     i__1 = k - 2;
00822                     for (j = 0; j <= i__1; ++j) {
00823                         slassq_(&j, &a[(j + 1) * lda], &c__1, &scale, &s);
00824 /*                    U at A(0,1) */
00825                     }
00826                     s += s;
00827 /*                 double s for the off diagonal elements */
00828                     i__1 = lda + 1;
00829                     slassq_(&k, a, &i__1, &scale, &s);
00830 /*                 tri L at A(0,0) */
00831                     i__1 = k - 1;
00832                     i__2 = lda + 1;
00833                     slassq_(&i__1, &a[lda], &i__2, &scale, &s);
00834 /*                 tri U at A(0,1) */
00835                 }
00836             } else {
00837 /*              A is xpose */
00838                 if (ilu == 0) {
00839 /*                 A' is upper */
00840                     i__1 = k - 2;
00841                     for (j = 1; j <= i__1; ++j) {
00842                         slassq_(&j, &a[(k + j) * lda], &c__1, &scale, &s);
00843 /*                    U at A(0,k) */
00844                     }
00845                     i__1 = k - 2;
00846                     for (j = 0; j <= i__1; ++j) {
00847                         slassq_(&k, &a[j * lda], &c__1, &scale, &s);
00848 /*                    k by k-1 rect. at A(0,0) */
00849                     }
00850                     i__1 = k - 2;
00851                     for (j = 0; j <= i__1; ++j) {
00852                         i__2 = k - j - 1;
00853                         slassq_(&i__2, &a[j + 1 + (j + k - 1) * lda], &c__1, &
00854                                 scale, &s);
00855 /*                    L at A(0,k-1) */
00856                     }
00857                     s += s;
00858 /*                 double s for the off diagonal elements */
00859                     i__1 = k - 1;
00860                     i__2 = lda + 1;
00861                     slassq_(&i__1, &a[k * lda], &i__2, &scale, &s);
00862 /*                 tri U at A(0,k) */
00863                     i__1 = lda + 1;
00864                     slassq_(&k, &a[(k - 1) * lda], &i__1, &scale, &s);
00865 /*                 tri L at A(0,k-1) */
00866                 } else {
00867 /*                 A' is lower */
00868                     i__1 = k - 1;
00869                     for (j = 1; j <= i__1; ++j) {
00870                         slassq_(&j, &a[j * lda], &c__1, &scale, &s);
00871 /*                    U at A(0,0) */
00872                     }
00873                     i__1 = *n - 1;
00874                     for (j = k; j <= i__1; ++j) {
00875                         slassq_(&k, &a[j * lda], &c__1, &scale, &s);
00876 /*                    k by k-1 rect. at A(0,k) */
00877                     }
00878                     i__1 = k - 3;
00879                     for (j = 0; j <= i__1; ++j) {
00880                         i__2 = k - j - 2;
00881                         slassq_(&i__2, &a[j + 2 + j * lda], &c__1, &scale, &s)
00882                                 ;
00883 /*                    L at A(1,0) */
00884                     }
00885                     s += s;
00886 /*                 double s for the off diagonal elements */
00887                     i__1 = lda + 1;
00888                     slassq_(&k, a, &i__1, &scale, &s);
00889 /*                 tri U at A(0,0) */
00890                     i__1 = k - 1;
00891                     i__2 = lda + 1;
00892                     slassq_(&i__1, &a[1], &i__2, &scale, &s);
00893 /*                 tri L at A(1,0) */
00894                 }
00895             }
00896         } else {
00897 /*           n is even */
00898             if (ifm == 1) {
00899 /*              A is normal */
00900                 if (ilu == 0) {
00901 /*                 A is upper */
00902                     i__1 = k - 2;
00903                     for (j = 0; j <= i__1; ++j) {
00904                         i__2 = k - j - 1;
00905                         slassq_(&i__2, &a[k + j + 2 + j * lda], &c__1, &scale, 
00906                                  &s);
00907 /*                    L at A(k+1,0) */
00908                     }
00909                     i__1 = k - 1;
00910                     for (j = 0; j <= i__1; ++j) {
00911                         i__2 = k + j;
00912                         slassq_(&i__2, &a[j * lda], &c__1, &scale, &s);
00913 /*                    trap U at A(0,0) */
00914                     }
00915                     s += s;
00916 /*                 double s for the off diagonal elements */
00917                     i__1 = lda + 1;
00918                     slassq_(&k, &a[k + 1], &i__1, &scale, &s);
00919 /*                 tri L at A(k+1,0) */
00920                     i__1 = lda + 1;
00921                     slassq_(&k, &a[k], &i__1, &scale, &s);
00922 /*                 tri U at A(k,0) */
00923                 } else {
00924 /*                 ilu=1 & A is lower */
00925                     i__1 = k - 1;
00926                     for (j = 0; j <= i__1; ++j) {
00927                         i__2 = *n - j - 1;
00928                         slassq_(&i__2, &a[j + 2 + j * lda], &c__1, &scale, &s)
00929                                 ;
00930 /*                    trap L at A(1,0) */
00931                     }
00932                     i__1 = k - 1;
00933                     for (j = 1; j <= i__1; ++j) {
00934                         slassq_(&j, &a[j * lda], &c__1, &scale, &s);
00935 /*                    U at A(0,0) */
00936                     }
00937                     s += s;
00938 /*                 double s for the off diagonal elements */
00939                     i__1 = lda + 1;
00940                     slassq_(&k, &a[1], &i__1, &scale, &s);
00941 /*                 tri L at A(1,0) */
00942                     i__1 = lda + 1;
00943                     slassq_(&k, a, &i__1, &scale, &s);
00944 /*                 tri U at A(0,0) */
00945                 }
00946             } else {
00947 /*              A is xpose */
00948                 if (ilu == 0) {
00949 /*                 A' is upper */
00950                     i__1 = k - 1;
00951                     for (j = 1; j <= i__1; ++j) {
00952                         slassq_(&j, &a[(k + 1 + j) * lda], &c__1, &scale, &s);
00953 /*                    U at A(0,k+1) */
00954                     }
00955                     i__1 = k - 1;
00956                     for (j = 0; j <= i__1; ++j) {
00957                         slassq_(&k, &a[j * lda], &c__1, &scale, &s);
00958 /*                    k by k rect. at A(0,0) */
00959                     }
00960                     i__1 = k - 2;
00961                     for (j = 0; j <= i__1; ++j) {
00962                         i__2 = k - j - 1;
00963                         slassq_(&i__2, &a[j + 1 + (j + k) * lda], &c__1, &
00964                                 scale, &s);
00965 /*                    L at A(0,k) */
00966                     }
00967                     s += s;
00968 /*                 double s for the off diagonal elements */
00969                     i__1 = lda + 1;
00970                     slassq_(&k, &a[(k + 1) * lda], &i__1, &scale, &s);
00971 /*                 tri U at A(0,k+1) */
00972                     i__1 = lda + 1;
00973                     slassq_(&k, &a[k * lda], &i__1, &scale, &s);
00974 /*                 tri L at A(0,k) */
00975                 } else {
00976 /*                 A' is lower */
00977                     i__1 = k - 1;
00978                     for (j = 1; j <= i__1; ++j) {
00979                         slassq_(&j, &a[(j + 1) * lda], &c__1, &scale, &s);
00980 /*                    U at A(0,1) */
00981                     }
00982                     i__1 = *n;
00983                     for (j = k + 1; j <= i__1; ++j) {
00984                         slassq_(&k, &a[j * lda], &c__1, &scale, &s);
00985 /*                    k by k rect. at A(0,k+1) */
00986                     }
00987                     i__1 = k - 2;
00988                     for (j = 0; j <= i__1; ++j) {
00989                         i__2 = k - j - 1;
00990                         slassq_(&i__2, &a[j + 1 + j * lda], &c__1, &scale, &s)
00991                                 ;
00992 /*                    L at A(0,0) */
00993                     }
00994                     s += s;
00995 /*                 double s for the off diagonal elements */
00996                     i__1 = lda + 1;
00997                     slassq_(&k, &a[lda], &i__1, &scale, &s);
00998 /*                 tri L at A(0,1) */
00999                     i__1 = lda + 1;
01000                     slassq_(&k, a, &i__1, &scale, &s);
01001 /*                 tri U at A(0,0) */
01002                 }
01003             }
01004         }
01005         value = scale * sqrt(s);
01006     }
01007 
01008     ret_val = value;
01009     return ret_val;
01010 
01011 /*     End of SLANSF */
01012 
01013 } /* slansf_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:10