00001 SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO )
00002
00003
00004
00005
00006
00007
00008
00009 CHARACTER UPLO
00010 INTEGER INFO, KB, LDA, LDW, N, NB
00011
00012
00013 INTEGER IPIV( * )
00014 DOUBLE PRECISION A( LDA, * ), W( LDW, * )
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098 DOUBLE PRECISION ZERO, ONE
00099 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00100 DOUBLE PRECISION EIGHT, SEVTEN
00101 PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
00102
00103
00104 INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
00105 $ KSTEP, KW
00106 DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1,
00107 $ ROWMAX, T
00108
00109
00110 LOGICAL LSAME
00111 INTEGER IDAMAX
00112 EXTERNAL LSAME, IDAMAX
00113
00114
00115 EXTERNAL DCOPY, DGEMM, DGEMV, DSCAL, DSWAP
00116
00117
00118 INTRINSIC ABS, MAX, MIN, SQRT
00119
00120
00121
00122 INFO = 0
00123
00124
00125
00126 ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
00127
00128 IF( LSAME( UPLO, 'U' ) ) THEN
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138 K = N
00139 10 CONTINUE
00140 KW = NB + K - N
00141
00142
00143
00144 IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
00145 $ GO TO 30
00146
00147
00148
00149 CALL DCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 )
00150 IF( K.LT.N )
00151 $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ), LDA,
00152 $ W( K, KW+1 ), LDW, ONE, W( 1, KW ), 1 )
00153
00154 KSTEP = 1
00155
00156
00157
00158
00159 ABSAKK = ABS( W( K, KW ) )
00160
00161
00162
00163
00164 IF( K.GT.1 ) THEN
00165 IMAX = IDAMAX( K-1, W( 1, KW ), 1 )
00166 COLMAX = ABS( W( IMAX, KW ) )
00167 ELSE
00168 COLMAX = ZERO
00169 END IF
00170
00171 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00172
00173
00174
00175 IF( INFO.EQ.0 )
00176 $ INFO = K
00177 KP = K
00178 ELSE
00179 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00180
00181
00182
00183 KP = K
00184 ELSE
00185
00186
00187
00188 CALL DCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
00189 CALL DCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
00190 $ W( IMAX+1, KW-1 ), 1 )
00191 IF( K.LT.N )
00192 $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ),
00193 $ LDA, W( IMAX, KW+1 ), LDW, ONE,
00194 $ W( 1, KW-1 ), 1 )
00195
00196
00197
00198
00199 JMAX = IMAX + IDAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 )
00200 ROWMAX = ABS( W( JMAX, KW-1 ) )
00201 IF( IMAX.GT.1 ) THEN
00202 JMAX = IDAMAX( IMAX-1, W( 1, KW-1 ), 1 )
00203 ROWMAX = MAX( ROWMAX, ABS( W( JMAX, KW-1 ) ) )
00204 END IF
00205
00206 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00207
00208
00209
00210 KP = K
00211 ELSE IF( ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN
00212
00213
00214
00215
00216 KP = IMAX
00217
00218
00219
00220 CALL DCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 )
00221 ELSE
00222
00223
00224
00225
00226 KP = IMAX
00227 KSTEP = 2
00228 END IF
00229 END IF
00230
00231 KK = K - KSTEP + 1
00232 KKW = NB + KK - N
00233
00234
00235
00236 IF( KP.NE.KK ) THEN
00237
00238
00239
00240 A( KP, K ) = A( KK, K )
00241 CALL DCOPY( K-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ),
00242 $ LDA )
00243 CALL DCOPY( KP, A( 1, KK ), 1, A( 1, KP ), 1 )
00244
00245
00246
00247 CALL DSWAP( N-KK+1, A( KK, KK ), LDA, A( KP, KK ), LDA )
00248 CALL DSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ),
00249 $ LDW )
00250 END IF
00251
00252 IF( KSTEP.EQ.1 ) THEN
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 CALL DCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 )
00263 R1 = ONE / A( K, K )
00264 CALL DSCAL( K-1, R1, A( 1, K ), 1 )
00265 ELSE
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275 IF( K.GT.2 ) THEN
00276
00277
00278
00279 D21 = W( K-1, KW )
00280 D11 = W( K, KW ) / D21
00281 D22 = W( K-1, KW-1 ) / D21
00282 T = ONE / ( D11*D22-ONE )
00283 D21 = T / D21
00284 DO 20 J = 1, K - 2
00285 A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) )
00286 A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) )
00287 20 CONTINUE
00288 END IF
00289
00290
00291
00292 A( K-1, K-1 ) = W( K-1, KW-1 )
00293 A( K-1, K ) = W( K-1, KW )
00294 A( K, K ) = W( K, KW )
00295 END IF
00296 END IF
00297
00298
00299
00300 IF( KSTEP.EQ.1 ) THEN
00301 IPIV( K ) = KP
00302 ELSE
00303 IPIV( K ) = -KP
00304 IPIV( K-1 ) = -KP
00305 END IF
00306
00307
00308
00309 K = K - KSTEP
00310 GO TO 10
00311
00312 30 CONTINUE
00313
00314
00315
00316
00317
00318
00319
00320 DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB
00321 JB = MIN( NB, K-J+1 )
00322
00323
00324
00325 DO 40 JJ = J, J + JB - 1
00326 CALL DGEMV( 'No transpose', JJ-J+1, N-K, -ONE,
00327 $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE,
00328 $ A( J, JJ ), 1 )
00329 40 CONTINUE
00330
00331
00332
00333 CALL DGEMM( 'No transpose', 'Transpose', J-1, JB, N-K, -ONE,
00334 $ A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, ONE,
00335 $ A( 1, J ), LDA )
00336 50 CONTINUE
00337
00338
00339
00340
00341 J = K + 1
00342 60 CONTINUE
00343 JJ = J
00344 JP = IPIV( J )
00345 IF( JP.LT.0 ) THEN
00346 JP = -JP
00347 J = J + 1
00348 END IF
00349 J = J + 1
00350 IF( JP.NE.JJ .AND. J.LE.N )
00351 $ CALL DSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA )
00352 IF( J.LE.N )
00353 $ GO TO 60
00354
00355
00356
00357 KB = N - K
00358
00359 ELSE
00360
00361
00362
00363
00364
00365
00366
00367 K = 1
00368 70 CONTINUE
00369
00370
00371
00372 IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N )
00373 $ GO TO 90
00374
00375
00376
00377 CALL DCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 )
00378 CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ), LDA,
00379 $ W( K, 1 ), LDW, ONE, W( K, K ), 1 )
00380
00381 KSTEP = 1
00382
00383
00384
00385
00386 ABSAKK = ABS( W( K, K ) )
00387
00388
00389
00390
00391 IF( K.LT.N ) THEN
00392 IMAX = K + IDAMAX( N-K, W( K+1, K ), 1 )
00393 COLMAX = ABS( W( IMAX, K ) )
00394 ELSE
00395 COLMAX = ZERO
00396 END IF
00397
00398 IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
00399
00400
00401
00402 IF( INFO.EQ.0 )
00403 $ INFO = K
00404 KP = K
00405 ELSE
00406 IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
00407
00408
00409
00410 KP = K
00411 ELSE
00412
00413
00414
00415 CALL DCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 )
00416 CALL DCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ),
00417 $ 1 )
00418 CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ),
00419 $ LDA, W( IMAX, 1 ), LDW, ONE, W( K, K+1 ), 1 )
00420
00421
00422
00423
00424 JMAX = K - 1 + IDAMAX( IMAX-K, W( K, K+1 ), 1 )
00425 ROWMAX = ABS( W( JMAX, K+1 ) )
00426 IF( IMAX.LT.N ) THEN
00427 JMAX = IMAX + IDAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 )
00428 ROWMAX = MAX( ROWMAX, ABS( W( JMAX, K+1 ) ) )
00429 END IF
00430
00431 IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
00432
00433
00434
00435 KP = K
00436 ELSE IF( ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN
00437
00438
00439
00440
00441 KP = IMAX
00442
00443
00444
00445 CALL DCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 )
00446 ELSE
00447
00448
00449
00450
00451 KP = IMAX
00452 KSTEP = 2
00453 END IF
00454 END IF
00455
00456 KK = K + KSTEP - 1
00457
00458
00459
00460 IF( KP.NE.KK ) THEN
00461
00462
00463
00464 A( KP, K ) = A( KK, K )
00465 CALL DCOPY( KP-K-1, A( K+1, KK ), 1, A( KP, K+1 ), LDA )
00466 CALL DCOPY( N-KP+1, A( KP, KK ), 1, A( KP, KP ), 1 )
00467
00468
00469
00470 CALL DSWAP( KK, A( KK, 1 ), LDA, A( KP, 1 ), LDA )
00471 CALL DSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW )
00472 END IF
00473
00474 IF( KSTEP.EQ.1 ) THEN
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484 CALL DCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 )
00485 IF( K.LT.N ) THEN
00486 R1 = ONE / A( K, K )
00487 CALL DSCAL( N-K, R1, A( K+1, K ), 1 )
00488 END IF
00489 ELSE
00490
00491
00492
00493
00494
00495
00496
00497
00498 IF( K.LT.N-1 ) THEN
00499
00500
00501
00502 D21 = W( K+1, K )
00503 D11 = W( K+1, K+1 ) / D21
00504 D22 = W( K, K ) / D21
00505 T = ONE / ( D11*D22-ONE )
00506 D21 = T / D21
00507 DO 80 J = K + 2, N
00508 A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) )
00509 A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) )
00510 80 CONTINUE
00511 END IF
00512
00513
00514
00515 A( K, K ) = W( K, K )
00516 A( K+1, K ) = W( K+1, K )
00517 A( K+1, K+1 ) = W( K+1, K+1 )
00518 END IF
00519 END IF
00520
00521
00522
00523 IF( KSTEP.EQ.1 ) THEN
00524 IPIV( K ) = KP
00525 ELSE
00526 IPIV( K ) = -KP
00527 IPIV( K+1 ) = -KP
00528 END IF
00529
00530
00531
00532 K = K + KSTEP
00533 GO TO 70
00534
00535 90 CONTINUE
00536
00537
00538
00539
00540
00541
00542
00543 DO 110 J = K, N, NB
00544 JB = MIN( NB, N-J+1 )
00545
00546
00547
00548 DO 100 JJ = J, J + JB - 1
00549 CALL DGEMV( 'No transpose', J+JB-JJ, K-1, -ONE,
00550 $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE,
00551 $ A( JJ, JJ ), 1 )
00552 100 CONTINUE
00553
00554
00555
00556 IF( J+JB.LE.N )
00557 $ CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
00558 $ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW,
00559 $ ONE, A( J+JB, J ), LDA )
00560 110 CONTINUE
00561
00562
00563
00564
00565 J = K - 1
00566 120 CONTINUE
00567 JJ = J
00568 JP = IPIV( J )
00569 IF( JP.LT.0 ) THEN
00570 JP = -JP
00571 J = J - 1
00572 END IF
00573 J = J - 1
00574 IF( JP.NE.JJ .AND. J.GE.1 )
00575 $ CALL DSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA )
00576 IF( J.GE.1 )
00577 $ GO TO 120
00578
00579
00580
00581 KB = K - 1
00582
00583 END IF
00584 RETURN
00585
00586
00587
00588 END