You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

dlansf.f 33 kB

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