Skip to content

Commit 9317095

Browse files
iyamazakiiyamazaki
committed
fix corner-case (N=1) for sytrs_aa (also fix typos in documentation)
1 parent eecc00b commit 9317095

File tree

6 files changed

+394
-249
lines changed

6 files changed

+394
-249
lines changed

‎SRC/chetrs_aa.f

Lines changed: 71 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -200,24 +200,31 @@ SUBROUTINE CHETRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
200200
*
201201
IF( UPPER ) THEN
202202
*
203-
* Solve A*X = B, where A = U*T*U**T.
203+
* Solve A*X = B, where A = U**T*T*U.
204204
*
205-
* P**T * B
205+
* 1) Forward substitution with U**T
206206
*
207-
K = 1
208-
DO WHILE ( K.LE.N )
209-
KP = IPIV( K )
210-
IF( KP.NE.K )
211-
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
212-
K = K + 1
213-
END DO
207+
IF( N.GT.1 ) THEN
208+
*
209+
* Pivot, P**T * B -> B
210+
*
211+
K = 1
212+
DO WHILE ( K.LE.N )
213+
KP = IPIV( K )
214+
IF( KP.NE.K )
215+
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
216+
K = K + 1
217+
END DO
214218
*
215-
* Compute (U \P**T * B) -> B [ (U \P**T * B) ]
219+
* Compute U**T \ B -> B [ (U**T \P**T * B) ]
220+
*
221+
CALL CTRSM( 'L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ),
222+
$ LDA, B( 2, 1 ), LDB)
223+
END IF
216224
*
217-
CALL CTRSM('L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
218-
$ B( 2, 1 ), LDB)
225+
* 2) Solve with triangular matrix T
219226
*
220-
* Compute T \ B -> B [ T \ (U \P**T * B) ]
227+
* Compute T \ B -> B [ T \ (U**T \P**T * B) ]
221228
*
222229
CALL CLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1)
223230
IF( N.GT.1 ) THEN
@@ -228,65 +235,82 @@ SUBROUTINE CHETRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
228235
CALL CGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
229236
$ INFO)
230237
*
231-
* Compute (U**T \ B) -> B [ U**T \ (T \ (U \P**T * B) ) ]
238+
* 3) Backward substitution with U
232239
*
233-
CALL CTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
234-
$ B(2, 1), LDB)
240+
IF( N.GT.1 ) THEN
241+
*
242+
* Compute U \ B -> B [ U \ (T \ (U**T \P**T * B) ) ]
235243
*
236-
* Pivot, P * B [ P * (U**T \ (T \ (U \P**T * B) )) ]
244+
CALL CTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ),
245+
$ LDA, B(2, 1), LDB)
237246
*
238-
K = N
239-
DO WHILE ( K.GE.1 )
240-
KP = IPIV( K )
241-
IF( KP.NE.K )
242-
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
243-
K = K - 1
244-
END DO
247+
* Pivot, P * B -> B [ P * (U \ (T \ (U**T \P**T * B) )) ]
248+
*
249+
K = N
250+
DO WHILE ( K.GE.1 )
251+
KP = IPIV( K )
252+
IF( KP.NE.K )
253+
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
254+
K = K - 1
255+
END DO
256+
END IF
245257
*
246258
ELSE
247259
*
248260
* Solve A*X = B, where A = L*T*L**T.
249261
*
250-
* Pivot, P**T * B
262+
* 1) Forward substitution with L
251263
*
252-
K = 1
253-
DO WHILE ( K.LE.N )
254-
KP = IPIV( K )
255-
IF( KP.NE.K )
256-
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
257-
K = K + 1
258-
END DO
264+
IF( N.GT.1 ) THEN
265+
*
266+
* Pivot, P**T * B -> B
267+
*
268+
K = 1
269+
DO WHILE ( K.LE.N )
270+
KP = IPIV( K )
271+
IF( KP.NE.K )
272+
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
273+
K = K + 1
274+
END DO
259275
*
260-
* Compute (L \P**T * B) -> B [ (L \P**T * B) ]
276+
* Compute L \ B -> B [ (L \P**T * B) ]
277+
*
278+
CALL CTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1),
279+
$ LDA, B(2, 1), LDB )
280+
END IF
261281
*
262-
CALL CTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1), LDA,
263-
$ B(2, 1), LDB)
282+
* 2) Solve with triangular matrix T
264283
*
265284
* Compute T \ B -> B [ T \ (L \P**T * B) ]
266285
*
267286
CALL CLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1)
268287
IF( N.GT.1 ) THEN
269-
CALL CLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 1 ), 1)
288+
CALL CLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 1 ), 1 )
270289
CALL CLACPY( 'F', 1, N-1, A( 2, 1 ), LDA+1, WORK( 2*N ), 1)
271290
CALL CLACGV( N-1, WORK( 2*N ), 1 )
272291
END IF
273292
CALL CGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB,
274293
$ INFO)
275294
*
276-
* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ]
295+
* 3) Backward substitution with L**T
277296
*
278-
CALL CTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA,
279-
$ B( 2, 1 ), LDB)
297+
IF( N.GT.1 ) THEN
298+
*
299+
* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ]
280300
*
281-
* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ]
301+
CALL CTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ),
302+
$ LDA, B( 2, 1 ), LDB )
282303
*
283-
K = N
284-
DO WHILE ( K.GE.1 )
285-
KP = IPIV( K )
286-
IF( KP.NE.K )
287-
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
288-
K = K - 1
289-
END DO
304+
* Pivot, P * B -> B [ P * (L**T \ (T \ (L \P**T * B) )) ]
305+
*
306+
K = N
307+
DO WHILE ( K.GE.1 )
308+
KP = IPIV( K )
309+
IF( KP.NE.K )
310+
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
311+
K = K - 1
312+
END DO
313+
END IF
290314
*
291315
END IF
292316
*

‎SRC/csytrs_aa.f

Lines changed: 63 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -200,22 +200,29 @@ SUBROUTINE CSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
200200
*
201201
IF( UPPER ) THEN
202202
*
203-
* Solve A*X = B, where A = U*T*U**T.
203+
* Solve A*X = B, where A = U**T*T*U.
204204
*
205-
* Pivot, P**T * B
205+
* 1) Forward substitution with U**T
206206
*
207-
DO K = 1, N
208-
KP = IPIV( K )
209-
IF( KP.NE.K )
210-
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
211-
END DO
207+
IF( N.GT.1 ) THEN
208+
*
209+
* Pivot, P**T * B -> B
210+
*
211+
DO K = 1, N
212+
KP = IPIV( K )
213+
IF( KP.NE.K )
214+
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
215+
END DO
212216
*
213-
* Compute (U \P**T * B) -> B [ (U \P**T * B) ]
217+
* Compute U**T \ B -> B [ (U**T \P**T * B) ]
218+
*
219+
CALL CTRSM( 'L', 'U', 'T', 'U', N-1, NRHS, ONE, A( 1, 2 ),
220+
$ LDA, B( 2, 1 ), LDB)
221+
END IF
214222
*
215-
CALL CTRSM('L', 'U', 'T', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
216-
$ B( 2, 1 ), LDB)
223+
* 2) Solve with triangular matrix T
217224
*
218-
* Compute T \ B -> B [ T \ (U \P**T * B) ]
225+
* Compute T \ B -> B [ T \ (U**T \P**T * B) ]
219226
*
220227
CALL CLACPY( 'F', 1, N, A( 1, 1 ), LDA+1, WORK( N ), 1)
221228
IF( N.GT.1 ) THEN
@@ -225,35 +232,48 @@ SUBROUTINE CSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
225232
CALL CGTSV( N, NRHS, WORK( 1 ), WORK( N ), WORK( 2*N ), B, LDB,
226233
$ INFO )
227234
*
228-
* Compute (U**T \ B) -> B [ U**T \ (T \ (U \P**T * B) ) ]
235+
* 3) Backward substitution with U
229236
*
230-
CALL CTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA,
231-
$ B( 2, 1 ), LDB)
237+
IF( N.GT.1 ) THEN
238+
*
239+
* Compute U \ B -> B [ U \ (T \ (U**T \P**T * B) ) ]
232240
*
233-
* Pivot, P * B [ P * (U**T \ (T \ (U \P**T * B) )) ]
241+
CALL CTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ),
242+
$ LDA, B( 2, 1 ), LDB)
234243
*
235-
DO K = N, 1, -1
236-
KP = IPIV( K )
237-
IF( KP.NE.K )
238-
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
239-
END DO
244+
* Pivot, P * B -> B [ P * (U**T \ (T \ (U \P**T * B) )) ]
245+
*
246+
DO K = N, 1, -1
247+
KP = IPIV( K )
248+
IF( KP.NE.K )
249+
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
250+
END DO
251+
END IF
240252
*
241253
ELSE
242254
*
243255
* Solve A*X = B, where A = L*T*L**T.
244256
*
245-
* Pivot, P**T * B
257+
* 1) Forward substitution with L
246258
*
247-
DO K = 1, N
248-
KP = IPIV( K )
249-
IF( KP.NE.K )
250-
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
251-
END DO
259+
IF( N.GT.1 ) THEN
252260
*
253-
* Compute (L \P**T * B) -> B [ (L \P**T * B) ]
261+
* Pivot, P**T * B -> B
262+
*
263+
DO K = 1, N
264+
KP = IPIV( K )
265+
IF( KP.NE.K )
266+
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
267+
END DO
268+
*
269+
* Compute L \ B -> B [ (L \P**T * B) ]
270+
*
271+
CALL CTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ),
272+
$ LDA, B( 2, 1 ), LDB)
273+
END IF
274+
*
275+
* 2) Solve with triangular matrix T
254276
*
255-
CALL CTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA,
256-
$ B( 2, 1 ), LDB)
257277
*
258278
* Compute T \ B -> B [ T \ (L \P**T * B) ]
259279
*
@@ -265,18 +285,23 @@ SUBROUTINE CSYTRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB,
265285
CALL CGTSV( N, NRHS, WORK( 1 ), WORK(N), WORK( 2*N ), B, LDB,
266286
$ INFO)
267287
*
268-
* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ]
288+
* 3) Backward substitution with L**T
269289
*
270-
CALL CTRSM( 'L', 'L', 'T', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA,
271-
$ B( 2, 1 ), LDB)
290+
IF( N.GT.1 ) THEN
291+
*
292+
* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ]
272293
*
273-
* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ]
294+
CALL CTRSM( 'L', 'L', 'T', 'U', N-1, NRHS, ONE, A( 2, 1 ),
295+
$ LDA, B( 2, 1 ), LDB)
274296
*
275-
DO K = N, 1, -1
276-
KP = IPIV( K )
277-
IF( KP.NE.K )
278-
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
279-
END DO
297+
* Pivot, P * B -> B [ P * (L**T \ (T \ (L \P**T * B) )) ]
298+
*
299+
DO K = N, 1, -1
300+
KP = IPIV( K )
301+
IF( KP.NE.K )
302+
$ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
303+
END DO
304+
END IF
280305
*
281306
END IF
282307
*

0 commit comments

Comments
 (0)