LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ ctrsen()

subroutine ctrsen ( character  job,
character  compq,
logical, dimension( * )  select,
integer  n,
complex, dimension( ldt, * )  t,
integer  ldt,
complex, dimension( ldq, * )  q,
integer  ldq,
complex, dimension( * )  w,
integer  m,
real  s,
real  sep,
complex, dimension( * )  work,
integer  lwork,
integer  info 
)

CTRSEN

Download CTRSEN + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CTRSEN reorders the Schur factorization of a complex matrix
 A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in
 the leading positions on the diagonal of the upper triangular matrix
 T, and the leading columns of Q form an orthonormal basis of the
 corresponding right invariant subspace.

 Optionally the routine computes the reciprocal condition numbers of
 the cluster of eigenvalues and/or the invariant subspace.
Parameters
[in]JOB
          JOB is CHARACTER*1
          Specifies whether condition numbers are required for the
          cluster of eigenvalues (S) or the invariant subspace (SEP):
          = 'N': none;
          = 'E': for eigenvalues only (S);
          = 'V': for invariant subspace only (SEP);
          = 'B': for both eigenvalues and invariant subspace (S and
                 SEP).
[in]COMPQ
          COMPQ is CHARACTER*1
          = 'V': update the matrix Q of Schur vectors;
          = 'N': do not update Q.
[in]SELECT
          SELECT is LOGICAL array, dimension (N)
          SELECT specifies the eigenvalues in the selected cluster. To
          select the j-th eigenvalue, SELECT(j) must be set to .TRUE..
[in]N
          N is INTEGER
          The order of the matrix T. N >= 0.
[in,out]T
          T is COMPLEX array, dimension (LDT,N)
          On entry, the upper triangular matrix T.
          On exit, T is overwritten by the reordered matrix T, with the
          selected eigenvalues as the leading diagonal elements.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T. LDT >= max(1,N).
[in,out]Q
          Q is COMPLEX array, dimension (LDQ,N)
          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
          On exit, if COMPQ = 'V', Q has been postmultiplied by the
          unitary transformation matrix which reorders T; the leading M
          columns of Q form an orthonormal basis for the specified
          invariant subspace.
          If COMPQ = 'N', Q is not referenced.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.
          LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
[out]W
          W is COMPLEX array, dimension (N)
          The reordered eigenvalues of T, in the same order as they
          appear on the diagonal of T.
[out]M
          M is INTEGER
          The dimension of the specified invariant subspace.
          0 <= M <= N.
[out]S
          S is REAL
          If JOB = 'E' or 'B', S is a lower bound on the reciprocal
          condition number for the selected cluster of eigenvalues.
          S cannot underestimate the true reciprocal condition number
          by more than a factor of sqrt(N). If M = 0 or N, S = 1.
          If JOB = 'N' or 'V', S is not referenced.
[out]SEP
          SEP is REAL
          If JOB = 'V' or 'B', SEP is the estimated reciprocal
          condition number of the specified invariant subspace. If
          M = 0 or N, SEP = norm(T).
          If JOB = 'N' or 'E', SEP is not referenced.
[out]WORK
          WORK is COMPLEX array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          If JOB = 'N', LWORK >= 1;
          if JOB = 'E', LWORK = max(1,M*(N-M));
          if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)).

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  CTRSEN first collects the selected eigenvalues by computing a unitary
  transformation Z to move them to the top left corner of T. In other
  words, the selected eigenvalues are the eigenvalues of T11 in:

          Z**H * T * Z = ( T11 T12 ) n1
                         (  0  T22 ) n2
                            n1  n2

  where N = n1+n2. The first
  n1 columns of Z span the specified invariant subspace of T.

  If T has been obtained from the Schur factorization of a matrix
  A = Q*T*Q**H, then the reordered Schur factorization of A is given by
  A = (Q*Z)*(Z**H*T*Z)*(Q*Z)**H, and the first n1 columns of Q*Z span the
  corresponding invariant subspace of A.

  The reciprocal condition number of the average of the eigenvalues of
  T11 may be returned in S. S lies between 0 (very badly conditioned)
  and 1 (very well conditioned). It is computed as follows. First we
  compute R so that

                         P = ( I  R ) n1
                             ( 0  0 ) n2
                               n1 n2

  is the projector on the invariant subspace associated with T11.
  R is the solution of the Sylvester equation:

                        T11*R - R*T22 = T12.

  Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
  the two-norm of M. Then S is computed as the lower bound

                      (1 + F-norm(R)**2)**(-1/2)

  on the reciprocal of 2-norm(P), the true reciprocal condition number.
  S cannot underestimate 1 / 2-norm(P) by more than a factor of
  sqrt(N).

  An approximate error bound for the computed average of the
  eigenvalues of T11 is

                         EPS * norm(T) / S

  where EPS is the machine precision.

  The reciprocal condition number of the right invariant subspace
  spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
  SEP is defined as the separation of T11 and T22:

                     sep( T11, T22 ) = sigma-min( C )

  where sigma-min(C) is the smallest singular value of the
  n1*n2-by-n1*n2 matrix

     C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )

  I(m) is an m by m identity matrix, and kprod denotes the Kronecker
  product. We estimate sigma-min(C) by the reciprocal of an estimate of
  the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
  cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).

  When SEP is small, small changes in T can cause large changes in
  the invariant subspace. An approximate bound on the maximum angular
  error in the computed right invariant subspace is

                      EPS * norm(T) / SEP

Definition at line 262 of file ctrsen.f.

264*
265* -- LAPACK computational routine --
266* -- LAPACK is a software package provided by Univ. of Tennessee, --
267* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268*
269* .. Scalar Arguments ..
270 CHARACTER COMPQ, JOB
271 INTEGER INFO, LDQ, LDT, LWORK, M, N
272 REAL S, SEP
273* ..
274* .. Array Arguments ..
275 LOGICAL SELECT( * )
276 COMPLEX Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
277* ..
278*
279* =====================================================================
280*
281* .. Parameters ..
282 REAL ZERO, ONE
283 parameter( zero = 0.0e+0, one = 1.0e+0 )
284* ..
285* .. Local Scalars ..
286 LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP
287 INTEGER IERR, K, KASE, KS, LWMIN, N1, N2, NN
288 REAL EST, RNORM, SCALE
289* ..
290* .. Local Arrays ..
291 INTEGER ISAVE( 3 )
292 REAL RWORK( 1 )
293* ..
294* .. External Functions ..
295 LOGICAL LSAME
296 REAL CLANGE, SROUNDUP_LWORK
297 EXTERNAL lsame, clange, sroundup_lwork
298* ..
299* .. External Subroutines ..
300 EXTERNAL clacn2, clacpy, ctrexc, ctrsyl, xerbla
301* ..
302* .. Intrinsic Functions ..
303 INTRINSIC max, sqrt
304* ..
305* .. Executable Statements ..
306*
307* Decode and test the input parameters.
308*
309 wantbh = lsame( job, 'B' )
310 wants = lsame( job, 'E' ) .OR. wantbh
311 wantsp = lsame( job, 'V' ) .OR. wantbh
312 wantq = lsame( compq, 'V' )
313*
314* Set M to the number of selected eigenvalues.
315*
316 m = 0
317 DO 10 k = 1, n
318 IF( SELECT( k ) )
319 $ m = m + 1
320 10 CONTINUE
321*
322 n1 = m
323 n2 = n - m
324 nn = n1*n2
325*
326 info = 0
327 lquery = ( lwork.EQ.-1 )
328*
329 IF( wantsp ) THEN
330 lwmin = max( 1, 2*nn )
331 ELSE IF( lsame( job, 'N' ) ) THEN
332 lwmin = 1
333 ELSE IF( lsame( job, 'E' ) ) THEN
334 lwmin = max( 1, nn )
335 END IF
336*
337 IF( .NOT.lsame( job, 'N' ) .AND. .NOT.wants .AND. .NOT.wantsp )
338 $ THEN
339 info = -1
340 ELSE IF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
341 info = -2
342 ELSE IF( n.LT.0 ) THEN
343 info = -4
344 ELSE IF( ldt.LT.max( 1, n ) ) THEN
345 info = -6
346 ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
347 info = -8
348 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
349 info = -14
350 END IF
351*
352 IF( info.EQ.0 ) THEN
353 work( 1 ) = sroundup_lwork(lwmin)
354 END IF
355*
356 IF( info.NE.0 ) THEN
357 CALL xerbla( 'CTRSEN', -info )
358 RETURN
359 ELSE IF( lquery ) THEN
360 RETURN
361 END IF
362*
363* Quick return if possible
364*
365 IF( m.EQ.n .OR. m.EQ.0 ) THEN
366 IF( wants )
367 $ s = one
368 IF( wantsp )
369 $ sep = clange( '1', n, n, t, ldt, rwork )
370 GO TO 40
371 END IF
372*
373* Collect the selected eigenvalues at the top left corner of T.
374*
375 ks = 0
376 DO 20 k = 1, n
377 IF( SELECT( k ) ) THEN
378 ks = ks + 1
379*
380* Swap the K-th eigenvalue to position KS.
381*
382 IF( k.NE.ks )
383 $ CALL ctrexc( compq, n, t, ldt, q, ldq, k, ks, ierr )
384 END IF
385 20 CONTINUE
386*
387 IF( wants ) THEN
388*
389* Solve the Sylvester equation for R:
390*
391* T11*R - R*T22 = scale*T12
392*
393 CALL clacpy( 'F', n1, n2, t( 1, n1+1 ), ldt, work, n1 )
394 CALL ctrsyl( 'N', 'N', -1, n1, n2, t, ldt, t( n1+1, n1+1 ),
395 $ ldt, work, n1, scale, ierr )
396*
397* Estimate the reciprocal of the condition number of the cluster
398* of eigenvalues.
399*
400 rnorm = clange( 'F', n1, n2, work, n1, rwork )
401 IF( rnorm.EQ.zero ) THEN
402 s = one
403 ELSE
404 s = scale / ( sqrt( scale*scale / rnorm+rnorm )*
405 $ sqrt( rnorm ) )
406 END IF
407 END IF
408*
409 IF( wantsp ) THEN
410*
411* Estimate sep(T11,T22).
412*
413 est = zero
414 kase = 0
415 30 CONTINUE
416 CALL clacn2( nn, work( nn+1 ), work, est, kase, isave )
417 IF( kase.NE.0 ) THEN
418 IF( kase.EQ.1 ) THEN
419*
420* Solve T11*R - R*T22 = scale*X.
421*
422 CALL ctrsyl( 'N', 'N', -1, n1, n2, t, ldt,
423 $ t( n1+1, n1+1 ), ldt, work, n1, scale,
424 $ ierr )
425 ELSE
426*
427* Solve T11**H*R - R*T22**H = scale*X.
428*
429 CALL ctrsyl( 'C', 'C', -1, n1, n2, t, ldt,
430 $ t( n1+1, n1+1 ), ldt, work, n1, scale,
431 $ ierr )
432 END IF
433 GO TO 30
434 END IF
435*
436 sep = scale / est
437 END IF
438*
439 40 CONTINUE
440*
441* Copy reordered eigenvalues to W.
442*
443 DO 50 k = 1, n
444 w( k ) = t( k, k )
445 50 CONTINUE
446*
447 work( 1 ) = sroundup_lwork(lwmin)
448*
449 RETURN
450*
451* End of CTRSEN
452*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine clacn2(n, v, x, est, kase, isave)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition clacn2.f:133
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:115
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine ctrexc(compq, n, t, ldt, q, ldq, ifst, ilst, info)
CTREXC
Definition ctrexc.f:126
subroutine ctrsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
CTRSYL
Definition ctrsyl.f:157
Here is the call graph for this function:
Here is the caller graph for this function: