Home
Browse Files
  1       DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
  2 *
  3 *  -- LAPACK routine (version 3.3.1)                                    --
  4 *
  5 *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --
  6 *  -- April 2011                                                      --
  7 *
  8 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  9 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 10 *
 11 *     .. Scalar Arguments ..
 12       CHARACTER          NORM, TRANSR, UPLO
 13       INTEGER            N
 14 *     ..
 15 *     .. Array Arguments ..
 16       DOUBLE PRECISION   A( 0: * ), WORK( 0: * )
 17 *     ..
 18 *
 19 *  Purpose
 20 *  =======
 21 *
 22 *  DLANSF returns the value of the one norm, or the Frobenius norm, or
 23 *  the infinity norm, or the element of largest absolute value of a
 24 *  real symmetric matrix A in RFP format.
 25 *
 26 *  Description
 27 *  ===========
 28 *
 29 *  DLANSF returns the value
 30 *
 31 *     DLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm'
 32 *              (
 33 *              ( norm1(A),         NORM = '1', 'O' or 'o'
 34 *              (
 35 *              ( normI(A),         NORM = 'I' or 'i'
 36 *              (
 37 *              ( normF(A),         NORM = 'F', 'f', 'E' or 'e'
 38 *
 39 *  where  norm1  denotes the  one norm of a matrix (maximum column sum),
 40 *  normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
 41 *  normF  denotes the  Frobenius norm of a matrix (square root of sum of
 42 *  squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm.
 43 *
 44 *  Arguments
 45 *  =========
 46 *
 47 *  NORM    (input) CHARACTER*1
 48 *          Specifies the value to be returned in DLANSF as described
 49 *          above.
 50 *
 51 *  TRANSR  (input) CHARACTER*1
 52 *          Specifies whether the RFP format of A is normal or
 53 *          transposed format.
 54 *          = 'N':  RFP format is Normal;
 55 *          = 'T':  RFP format is Transpose.
 56 *
 57 *  UPLO    (input) CHARACTER*1
 58 *           On entry, UPLO specifies whether the RFP matrix A came from
 59 *           an upper or lower triangular matrix as follows:
 60 *           = 'U': RFP A came from an upper triangular matrix;
 61 *           = 'L': RFP A came from a lower triangular matrix.
 62 *
 63 *  N       (input) INTEGER
 64 *          The order of the matrix A. N >= 0. When N = 0, DLANSF is
 65 *          set to zero.
 66 *
 67 *  A       (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 );
 68 *          On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')
 69 *          part of the symmetric matrix A stored in RFP format. See the
 70 *          "Notes" below for more details.
 71 *          Unchanged on exit.
 72 *
 73 *  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
 74 *          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
 75 *          WORK is not referenced.
 76 *
 77 *  Further Details
 78 *  ===============
 79 *
 80 *  We first consider Rectangular Full Packed (RFP) Format when N is
 81 *  even. We give an example where N = 6.
 82 *
 83 *      AP is Upper             AP is Lower
 84 *
 85 *   00 01 02 03 04 05       00
 86 *      11 12 13 14 15       10 11
 87 *         22 23 24 25       20 21 22
 88 *            33 34 35       30 31 32 33
 89 *               44 45       40 41 42 43 44
 90 *                  55       50 51 52 53 54 55
 91 *
 92 *
 93 *  Let TRANSR = 'N'. RFP holds AP as follows:
 94 *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
 95 *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
 96 *  the transpose of the first three columns of AP upper.
 97 *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
 98 *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
 99 *  the transpose of the last three columns of AP lower.
100 *  This covers the case N even and TRANSR = 'N'.
101 *
102 *         RFP A                   RFP A
103 *
104 *        03 04 05                33 43 53
105 *        13 14 15                00 44 54
106 *        23 24 25                10 11 55
107 *        33 34 35                20 21 22
108 *        00 44 45                30 31 32
109 *        01 11 55                40 41 42
110 *        02 12 22                50 51 52
111 *
112 *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
113 *  transpose of RFP A above. One therefore gets:
114 *
115 *
116 *           RFP A                   RFP A
117 *
118 *     03 13 23 33 00 01 02    33 00 10 20 30 40 50
119 *     04 14 24 34 44 11 12    43 44 11 21 31 41 51
120 *     05 15 25 35 45 55 22    53 54 55 22 32 42 52
121 *
122 *
123 *  We then consider Rectangular Full Packed (RFP) Format when N is
124 *  odd. We give an example where N = 5.
125 *
126 *     AP is Upper                 AP is Lower
127 *
128 *   00 01 02 03 04              00
129 *      11 12 13 14              10 11
130 *         22 23 24              20 21 22
131 *            33 34              30 31 32 33
132 *               44              40 41 42 43 44
133 *
134 *
135 *  Let TRANSR = 'N'. RFP holds AP as follows:
136 *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
137 *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
138 *  the transpose of the first two columns of AP upper.
139 *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
140 *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
141 *  the transpose of the last two columns of AP lower.
142 *  This covers the case N odd and TRANSR = 'N'.
143 *
144 *         RFP A                   RFP A
145 *
146 *        02 03 04                00 33 43
147 *        12 13 14                10 11 44
148 *        22 23 24                20 21 22
149 *        00 33 34                30 31 32
150 *        01 11 44                40 41 42
151 *
152 *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
153 *  transpose of RFP A above. One therefore gets:
154 *
155 *           RFP A                   RFP A
156 *
157 *     02 12 22 00 01             00 10 20 30 40 50
158 *     03 13 23 33 11             33 11 21 31 41 51
159 *     04 14 24 34 44             43 44 22 32 42 52
160 *
161 *  Reference
162 *  =========
163 *
164 *  =====================================================================
165 *
166 *     .. Parameters ..
167       DOUBLE PRECISION   ONE, ZERO
168       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
169 *     ..
170 *     .. Local Scalars ..
171       INTEGER            I, J, IFM, ILU, NOE, N1, K, L, LDA
172       DOUBLE PRECISION   SCALE, S, VALUE, AA
173 *     ..
174 *     .. External Functions ..
175       LOGICAL            LSAME
176       INTEGER            IDAMAX
177       EXTERNAL           LSAME, IDAMAX
178 *     ..
179 *     .. External Subroutines ..
180       EXTERNAL           DLASSQ
181 *     ..
182 *     .. Intrinsic Functions ..
183       INTRINSIC          ABS, MAX, SQRT
184 *     ..
185 *     .. Executable Statements ..
186 *
187       IF( N.EQ.0 ) THEN
188          DLANSF = ZERO
189          RETURN
190       END IF
191 *
192 *     set noe = 1 if n is odd. if n is even set noe=0
193 *
194       NOE = 1
195       IF( MOD( N, 2 ).EQ.0 )
196      $   NOE = 0
197 *
198 *     set ifm = 0 when form='T or 't' and 1 otherwise
199 *
200       IFM = 1
201       IF( LSAME( TRANSR, 'T' ) )
202      $   IFM = 0
203 *
204 *     set ilu = 0 when uplo='U or 'u' and 1 otherwise
205 *
206       ILU = 1
207       IF( LSAME( UPLO, 'U' ) )
208      $   ILU = 0
209 *
210 *     set lda = (n+1)/2 when ifm = 0
211 *     set lda = n when ifm = 1 and noe = 1
212 *     set lda = n+1 when ifm = 1 and noe = 0
213 *
214       IF( IFM.EQ.1 ) THEN
215          IF( NOE.EQ.1 ) THEN
216             LDA = N
217          ELSE
218 *           noe=0
219             LDA = N + 1
220          END IF
221       ELSE
222 *        ifm=0
223          LDA = ( N+1 ) / 2
224       END IF
225 *
226       IF( LSAME( NORM, 'M' ) ) THEN
227 *
228 *       Find max(abs(A(i,j))).
229 *
230          K = ( N+1 ) / 2
231          VALUE = ZERO
232          IF( NOE.EQ.1 ) THEN
233 *           n is odd
234             IF( IFM.EQ.1 ) THEN
235 *           A is n by k
236                DO J = 0, K - 1
237                   DO I = 0, N - 1
238                      VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
239                   END DO
240                END DO
241             ELSE
242 *              xpose case; A is k by n
243                DO J = 0, N - 1
244                   DO I = 0, K - 1
245                      VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
246                   END DO
247                END DO
248             END IF
249          ELSE
250 *           n is even
251             IF( IFM.EQ.1 ) THEN
252 *              A is n+1 by k
253                DO J = 0, K - 1
254                   DO I = 0, N
255                      VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
256                   END DO
257                END DO
258             ELSE
259 *              xpose case; A is k by n+1
260                DO J = 0, N
261                   DO I = 0, K - 1
262                      VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
263                   END DO
264                END DO
265             END IF
266          END IF
267       ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
268      $         ( NORM.EQ.'1' ) ) THEN
269 *
270 *        Find normI(A) ( = norm1(A), since A is symmetric).
271 *
272          IF( IFM.EQ.1 ) THEN
273             K = N / 2
274             IF( NOE.EQ.1 ) THEN
275 *              n is odd
276                IF( ILU.EQ.0 ) THEN
277                   DO I = 0, K - 1
278                      WORK( I ) = ZERO
279                   END DO
280                   DO J = 0, K
281                      S = ZERO
282                      DO I = 0, K + J - 1
283                         AA = ABS( A( I+J*LDA ) )
284 *                       -> A(i,j+k)
285                         S = S + AA
286                         WORK( I ) = WORK( I ) + AA
287                      END DO
288                      AA = ABS( A( I+J*LDA ) )
289 *                    -> A(j+k,j+k)
290                      WORK( J+K ) = S + AA
291                      IF( I.EQ.K+K )
292      $                  GO TO 10
293                      I = I + 1
294                      AA = ABS( A( I+J*LDA ) )
295 *                    -> A(j,j)
296                      WORK( J ) = WORK( J ) + AA
297                      S = ZERO
298                      DO L = J + 1, K - 1
299                         I = I + 1
300                         AA = ABS( A( I+J*LDA ) )
301 *                       -> A(l,j)
302                         S = S + AA
303                         WORK( L ) = WORK( L ) + AA
304                      END DO
305                      WORK( J ) = WORK( J ) + S
306                   END DO
307    10             CONTINUE
308                   I = IDAMAX( N, WORK, 1 )
309                   VALUE = WORK( I-1 )
310                ELSE
311 *                 ilu = 1
312                   K = K + 1
313 *                 k=(n+1)/2 for n odd and ilu=1
314                   DO I = K, N - 1
315                      WORK( I ) = ZERO
316                   END DO
317                   DO J = K - 1, 0, -1
318                      S = ZERO
319                      DO I = 0, J - 2
320                         AA = ABS( A( I+J*LDA ) )
321 *                       -> A(j+k,i+k)
322                         S = S + AA
323                         WORK( I+K ) = WORK( I+K ) + AA
324                      END DO
325                      IF( J.GT.0 ) THEN
326                         AA = ABS( A( I+J*LDA ) )
327 *                       -> A(j+k,j+k)
328                         S = S + AA
329                         WORK( I+K ) = WORK( I+K ) + S
330 *                       i=j
331                         I = I + 1
332                      END IF
333                      AA = ABS( A( I+J*LDA ) )
334 *                    -> A(j,j)
335                      WORK( J ) = AA
336                      S = ZERO
337                      DO L = J + 1, N - 1
338                         I = I + 1
339                         AA = ABS( A( I+J*LDA ) )
340 *                       -> A(l,j)
341                         S = S + AA
342                         WORK( L ) = WORK( L ) + AA
343                      END DO
344                      WORK( J ) = WORK( J ) + S
345                   END DO
346                   I = IDAMAX( N, WORK, 1 )
347                   VALUE = WORK( I-1 )
348                END IF
349             ELSE
350 *              n is even
351                IF( ILU.EQ.0 ) THEN
352                   DO I = 0, K - 1
353                      WORK( I ) = ZERO
354                   END DO
355                   DO J = 0, K - 1
356                      S = ZERO
357                      DO I = 0, K + J - 1
358                         AA = ABS( A( I+J*LDA ) )
359 *                       -> A(i,j+k)
360                         S = S + AA
361                         WORK( I ) = WORK( I ) + AA
362                      END DO
363                      AA = ABS( A( I+J*LDA ) )
364 *                    -> A(j+k,j+k)
365                      WORK( J+K ) = S + AA
366                      I = I + 1
367                      AA = ABS( A( I+J*LDA ) )
368 *                    -> A(j,j)
369                      WORK( J ) = WORK( J ) + AA
370                      S = ZERO
371                      DO L = J + 1, K - 1
372                         I = I + 1
373                         AA = ABS( A( I+J*LDA ) )
374 *                       -> A(l,j)
375                         S = S + AA
376                         WORK( L ) = WORK( L ) + AA
377                      END DO
378                      WORK( J ) = WORK( J ) + S
379                   END DO
380                   I = IDAMAX( N, WORK, 1 )
381                   VALUE = WORK( I-1 )
382                ELSE
383 *                 ilu = 1
384                   DO I = K, N - 1
385                      WORK( I ) = ZERO
386                   END DO
387                   DO J = K - 1, 0, -1
388                      S = ZERO
389                      DO I = 0, J - 1
390                         AA = ABS( A( I+J*LDA ) )
391 *                       -> A(j+k,i+k)
392                         S = S + AA
393                         WORK( I+K ) = WORK( I+K ) + AA
394                      END DO
395                      AA = ABS( A( I+J*LDA ) )
396 *                    -> A(j+k,j+k)
397                      S = S + AA
398                      WORK( I+K ) = WORK( I+K ) + S
399 *                    i=j
400                      I = I + 1
401                      AA = ABS( A( I+J*LDA ) )
402 *                    -> A(j,j)
403                      WORK( J ) = AA
404                      S = ZERO
405                      DO L = J + 1, N - 1
406                         I = I + 1
407                         AA = ABS( A( I+J*LDA ) )
408 *                       -> A(l,j)
409                         S = S + AA
410                         WORK( L ) = WORK( L ) + AA
411                      END DO
412                      WORK( J ) = WORK( J ) + S
413                   END DO
414                   I = IDAMAX( N, WORK, 1 )
415                   VALUE = WORK( I-1 )
416                END IF
417             END IF
418          ELSE
419 *           ifm=0
420             K = N / 2
421             IF( NOE.EQ.1 ) THEN
422 *              n is odd
423                IF( ILU.EQ.0 ) THEN
424                   N1 = K
425 *                 n/2
426                   K = K + 1
427 *                 k is the row size and lda
428                   DO I = N1, N - 1
429                      WORK( I ) = ZERO
430                   END DO
431                   DO J = 0, N1 - 1
432                      S = ZERO
433                      DO I = 0, K - 1
434                         AA = ABS( A( I+J*LDA ) )
435 *                       A(j,n1+i)
436                         WORK( I+N1 ) = WORK( I+N1 ) + AA
437                         S = S + AA
438                      END DO
439                      WORK( J ) = S
440                   END DO
441 *                 j=n1=k-1 is special
442                   S = ABS( A( 0+J*LDA ) )
443 *                 A(k-1,k-1)
444                   DO I = 1, K - 1
445                      AA = ABS( A( I+J*LDA ) )
446 *                    A(k-1,i+n1)
447                      WORK( I+N1 ) = WORK( I+N1 ) + AA
448                      S = S + AA
449                   END DO
450                   WORK( J ) = WORK( J ) + S
451                   DO J = K, N - 1
452                      S = ZERO
453                      DO I = 0, J - K - 1
454                         AA = ABS( A( I+J*LDA ) )
455 *                       A(i,j-k)
456                         WORK( I ) = WORK( I ) + AA
457                         S = S + AA
458                      END DO
459 *                    i=j-k
460                      AA = ABS( A( I+J*LDA ) )
461 *                    A(j-k,j-k)
462                      S = S + AA
463                      WORK( J-K ) = WORK( J-K ) + S
464                      I = I + 1
465                      S = ABS( A( I+J*LDA ) )
466 *                    A(j,j)
467                      DO L = J + 1, N - 1
468                         I = I + 1
469                         AA = ABS( A( I+J*LDA ) )
470 *                       A(j,l)
471                         WORK( L ) = WORK( L ) + AA
472                         S = S + AA
473                      END DO
474                      WORK( J ) = WORK( J ) + S
475                   END DO
476                   I = IDAMAX( N, WORK, 1 )
477                   VALUE = WORK( I-1 )
478                ELSE
479 *                 ilu=1
480                   K = K + 1
481 *                 k=(n+1)/2 for n odd and ilu=1
482                   DO I = K, N - 1
483                      WORK( I ) = ZERO
484                   END DO
485                   DO J = 0, K - 2
486 *                    process
487                      S = ZERO
488                      DO I = 0, J - 1
489                         AA = ABS( A( I+J*LDA ) )
490 *                       A(j,i)
491                         WORK( I ) = WORK( I ) + AA
492                         S = S + AA
493                      END DO
494                      AA = ABS( A( I+J*LDA ) )
495 *                    i=j so process of A(j,j)
496                      S = S + AA
497                      WORK( J ) = S
498 *                    is initialised here
499                      I = I + 1
500 *                    i=j process A(j+k,j+k)
501                      AA = ABS( A( I+J*LDA ) )
502                      S = AA
503                      DO L = K + J + 1, N - 1
504                         I = I + 1
505                         AA = ABS( A( I+J*LDA ) )
506 *                       A(l,k+j)
507                         S = S + AA
508                         WORK( L ) = WORK( L ) + AA
509                      END DO
510                      WORK( K+J ) = WORK( K+J ) + S
511                   END DO
512 *                 j=k-1 is special :process col A(k-1,0:k-1)
513                   S = ZERO
514                   DO I = 0, K - 2
515                      AA = ABS( A( I+J*LDA ) )
516 *                    A(k,i)
517                      WORK( I ) = WORK( I ) + AA
518                      S = S + AA
519                   END DO
520 *                 i=k-1
521                   AA = ABS( A( I+J*LDA ) )
522 *                 A(k-1,k-1)
523                   S = S + AA
524                   WORK( I ) = S
525 *                 done with col j=k+1
526                   DO J = K, N - 1
527 *                    process col j of A = A(j,0:k-1)
528                      S = ZERO
529                      DO I = 0, K - 1
530                         AA = ABS( A( I+J*LDA ) )
531 *                       A(j,i)
532                         WORK( I ) = WORK( I ) + AA
533                         S = S + AA
534                      END DO
535                      WORK( J ) = WORK( J ) + S
536                   END DO
537                   I = IDAMAX( N, WORK, 1 )
538                   VALUE = WORK( I-1 )
539                END IF
540             ELSE
541 *              n is even
542                IF( ILU.EQ.0 ) THEN
543                   DO I = K, N - 1
544                      WORK( I ) = ZERO
545                   END DO
546                   DO J = 0, K - 1
547                      S = ZERO
548                      DO I = 0, K - 1
549                         AA = ABS( A( I+J*LDA ) )
550 *                       A(j,i+k)
551                         WORK( I+K ) = WORK( I+K ) + AA
552                         S = S + AA
553                      END DO
554                      WORK( J ) = S
555                   END DO
556 *                 j=k
557                   AA = ABS( A( 0+J*LDA ) )
558 *                 A(k,k)
559                   S = AA
560                   DO I = 1, K - 1
561                      AA = ABS( A( I+J*LDA ) )
562 *                    A(k,k+i)
563                      WORK( I+K ) = WORK( I+K ) + AA
564                      S = S + AA
565                   END DO
566                   WORK( J ) = WORK( J ) + S
567                   DO J = K + 1, N - 1
568                      S = ZERO
569                      DO I = 0, J - 2 - K
570                         AA = ABS( A( I+J*LDA ) )
571 *                       A(i,j-k-1)
572                         WORK( I ) = WORK( I ) + AA
573                         S = S + AA
574                      END DO
575 *                     i=j-1-k
576                      AA = ABS( A( I+J*LDA ) )
577 *                    A(j-k-1,j-k-1)
578                      S = S + AA
579                      WORK( J-K-1 ) = WORK( J-K-1 ) + S
580                      I = I + 1
581                      AA = ABS( A( I+J*LDA ) )
582 *                    A(j,j)
583                      S = AA
584                      DO L = J + 1, N - 1
585                         I = I + 1
586                         AA = ABS( A( I+J*LDA ) )
587 *                       A(j,l)
588                         WORK( L ) = WORK( L ) + AA
589                         S = S + AA
590                      END DO
591                      WORK( J ) = WORK( J ) + S
592                   END DO
593 *                 j=n
594                   S = ZERO
595                   DO I = 0, K - 2
596                      AA = ABS( A( I+J*LDA ) )
597 *                    A(i,k-1)
598                      WORK( I ) = WORK( I ) + AA
599                      S = S + AA
600                   END DO
601 *                 i=k-1
602                   AA = ABS( A( I+J*LDA ) )
603 *                 A(k-1,k-1)
604                   S = S + AA
605                   WORK( I ) = WORK( I ) + S
606                   I = IDAMAX( N, WORK, 1 )
607                   VALUE = WORK( I-1 )
608                ELSE
609 *                 ilu=1
610                   DO I = K, N - 1
611                      WORK( I ) = ZERO
612                   END DO
613 *                 j=0 is special :process col A(k:n-1,k)
614                   S = ABS( A( 0 ) )
615 *                 A(k,k)
616                   DO I = 1, K - 1
617                      AA = ABS( A( I ) )
618 *                    A(k+i,k)
619                      WORK( I+K ) = WORK( I+K ) + AA
620                      S = S + AA
621                   END DO
622                   WORK( K ) = WORK( K ) + S
623                   DO J = 1, K - 1
624 *                    process
625                      S = ZERO
626                      DO I = 0, J - 2
627                         AA = ABS( A( I+J*LDA ) )
628 *                       A(j-1,i)
629                         WORK( I ) = WORK( I ) + AA
630                         S = S + AA
631                      END DO
632                      AA = ABS( A( I+J*LDA ) )
633 *                    i=j-1 so process of A(j-1,j-1)
634                      S = S + AA
635                      WORK( J-1 ) = S
636 *                    is initialised here
637                      I = I + 1
638 *                    i=j process A(j+k,j+k)
639                      AA = ABS( A( I+J*LDA ) )
640                      S = AA
641                      DO L = K + J + 1, N - 1
642                         I = I + 1
643                         AA = ABS( A( I+J*LDA ) )
644 *                       A(l,k+j)
645                         S = S + AA
646                         WORK( L ) = WORK( L ) + AA
647                      END DO
648                      WORK( K+J ) = WORK( K+J ) + S
649                   END DO
650 *                 j=k is special :process col A(k,0:k-1)
651                   S = ZERO
652                   DO I = 0, K - 2
653                      AA = ABS( A( I+J*LDA ) )
654 *                    A(k,i)
655                      WORK( I ) = WORK( I ) + AA
656                      S = S + AA
657                   END DO
658 *                 i=k-1
659                   AA = ABS( A( I+J*LDA ) )
660 *                 A(k-1,k-1)
661                   S = S + AA
662                   WORK( I ) = S
663 *                 done with col j=k+1
664                   DO J = K + 1, N
665 *                    process col j-1 of A = A(j-1,0:k-1)
666                      S = ZERO
667                      DO I = 0, K - 1
668                         AA = ABS( A( I+J*LDA ) )
669 *                       A(j-1,i)
670                         WORK( I ) = WORK( I ) + AA
671                         S = S + AA
672                      END DO
673                      WORK( J-1 ) = WORK( J-1 ) + S
674                   END DO
675                   I = IDAMAX( N, WORK, 1 )
676                   VALUE = WORK( I-1 )
677                END IF
678             END IF
679          END IF
680       ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
681 *
682 *       Find normF(A).
683 *
684          K = ( N+1 ) / 2
685          SCALE = ZERO
686          S = ONE
687          IF( NOE.EQ.1 ) THEN
688 *           n is odd
689             IF( IFM.EQ.1 ) THEN
690 *              A is normal
691                IF( ILU.EQ.0 ) THEN
692 *                 A is upper
693                   DO J = 0, K - 3
694                      CALL DLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S )
695 *                    L at A(k,0)
696                   END DO
697                   DO J = 0, K - 1
698                      CALL DLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S )
699 *                    trap U at A(0,0)
700                   END DO
701                   S = S + S
702 *                 double s for the off diagonal elements
703                   CALL DLASSQ( K-1, A( K ), LDA+1, SCALE, S )
704 *                 tri L at A(k,0)
705                   CALL DLASSQ( K, A( K-1 ), LDA+1, SCALE, S )
706 *                 tri U at A(k-1,0)
707                ELSE
708 *                 ilu=1 & A is lower
709                   DO J = 0, K - 1
710                      CALL DLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S )
711 *                    trap L at A(0,0)
712                   END DO
713                   DO J = 0, K - 2
714                      CALL DLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S )
715 *                    U at A(0,1)
716                   END DO
717                   S = S + S
718 *                 double s for the off diagonal elements
719                   CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
720 *                 tri L at A(0,0)
721                   CALL DLASSQ( K-1, A( 0+LDA ), LDA+1, SCALE, S )
722 *                 tri U at A(0,1)
723                END IF
724             ELSE
725 *              A is xpose
726                IF( ILU.EQ.0 ) THEN
727 *                 A**T is upper
728                   DO J = 1, K - 2
729                      CALL DLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S )
730 *                    U at A(0,k)
731                   END DO
732                   DO J = 0, K - 2
733                      CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
734 *                    k by k-1 rect. at A(0,0)
735                   END DO
736                   DO J = 0, K - 2
737                      CALL DLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1,
738      $                            SCALE, S )
739 *                    L at A(0,k-1)
740                   END DO
741                   S = S + S
742 *                 double s for the off diagonal elements
743                   CALL DLASSQ( K-1, A( 0+K*LDA ), LDA+1, SCALE, S )
744 *                 tri U at A(0,k)
745                   CALL DLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S )
746 *                 tri L at A(0,k-1)
747                ELSE
748 *                 A**T is lower
749                   DO J = 1, K - 1
750                      CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
751 *                    U at A(0,0)
752                   END DO
753                   DO J = K, N - 1
754                      CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
755 *                    k by k-1 rect. at A(0,k)
756                   END DO
757                   DO J = 0, K - 3
758                      CALL DLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S )
759 *                    L at A(1,0)
760                   END DO
761                   S = S + S
762 *                 double s for the off diagonal elements
763                   CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
764 *                 tri U at A(0,0)
765                   CALL DLASSQ( K-1, A( 1 ), LDA+1, SCALE, S )
766 *                 tri L at A(1,0)
767                END IF
768             END IF
769          ELSE
770 *           n is even
771             IF( IFM.EQ.1 ) THEN
772 *              A is normal
773                IF( ILU.EQ.0 ) THEN
774 *                 A is upper
775                   DO J = 0, K - 2
776                      CALL DLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S )
777 *                    L at A(k+1,0)
778                   END DO
779                   DO J = 0, K - 1
780                      CALL DLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S )
781 *                    trap U at A(0,0)
782                   END DO
783                   S = S + S
784 *                 double s for the off diagonal elements
785                   CALL DLASSQ( K, A( K+1 ), LDA+1, SCALE, S )
786 *                 tri L at A(k+1,0)
787                   CALL DLASSQ( K, A( K ), LDA+1, SCALE, S )
788 *                 tri U at A(k,0)
789                ELSE
790 *                 ilu=1 & A is lower
791                   DO J = 0, K - 1
792                      CALL DLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S )
793 *                    trap L at A(1,0)
794                   END DO
795                   DO J = 1, K - 1
796                      CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
797 *                    U at A(0,0)
798                   END DO
799                   S = S + S
800 *                 double s for the off diagonal elements
801                   CALL DLASSQ( K, A( 1 ), LDA+1, SCALE, S )
802 *                 tri L at A(1,0)
803                   CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
804 *                 tri U at A(0,0)
805                END IF
806             ELSE
807 *              A is xpose
808                IF( ILU.EQ.0 ) THEN
809 *                 A**T is upper
810                   DO J = 1, K - 1
811                      CALL DLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S )
812 *                    U at A(0,k+1)
813                   END DO
814                   DO J = 0, K - 1
815                      CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
816 *                    k by k rect. at A(0,0)
817                   END DO
818                   DO J = 0, K - 2
819                      CALL DLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE,
820      $                            S )
821 *                    L at A(0,k)
822                   END DO
823                   S = S + S
824 *                 double s for the off diagonal elements
825                   CALL DLASSQ( K, A( 0+( K+1 )*LDA ), LDA+1, SCALE, S )
826 *                 tri U at A(0,k+1)
827                   CALL DLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S )
828 *                 tri L at A(0,k)
829                ELSE
830 *                 A**T is lower
831                   DO J = 1, K - 1
832                      CALL DLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S )
833 *                    U at A(0,1)
834                   END DO
835                   DO J = K + 1, N
836                      CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
837 *                    k by k rect. at A(0,k+1)
838                   END DO
839                   DO J = 0, K - 2
840                      CALL DLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S )
841 *                    L at A(0,0)
842                   END DO
843                   S = S + S
844 *                 double s for the off diagonal elements
845                   CALL DLASSQ( K, A( LDA ), LDA+1, SCALE, S )
846 *                 tri L at A(0,1)
847                   CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
848 *                 tri U at A(0,0)
849                END IF
850             END IF
851          END IF
852          VALUE = SCALE*SQRT( S )
853       END IF
854 *
855       DLANSF = VALUE
856       RETURN
857 *
858 *     End of DLANSF
859 *
860       END
Home