 LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 subroutine ztgsyl ( character TRANS, integer IJOB, integer M, integer N, complex*16, dimension( lda, * ) A, integer LDA, complex*16, dimension( ldb, * ) B, integer LDB, complex*16, dimension( ldc, * ) C, integer LDC, complex*16, dimension( ldd, * ) D, integer LDD, complex*16, dimension( lde, * ) E, integer LDE, complex*16, dimension( ldf, * ) F, integer LDF, double precision SCALE, double precision DIF, complex*16, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer INFO )

ZTGSYL

Purpose:
ZTGSYL solves the generalized Sylvester equation:

A * R - L * B = scale * C            (1)
D * R - L * E = scale * F

where R and L are unknown m-by-n matrices, (A, D), (B, E) and
(C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
respectively, with complex entries. A, B, D and E are upper
triangular (i.e., (A,D) and (B,E) in generalized Schur form).

The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1
is an output scaling factor chosen to avoid overflow.

In matrix notation (1) is equivalent to solve Zx = scale*b, where Z
is defined as

Z = [ kron(In, A)  -kron(B**H, Im) ]        (2)
[ kron(In, D)  -kron(E**H, Im) ],

Here Ix is the identity matrix of size x and X**H is the conjugate
transpose of X. Kron(X, Y) is the Kronecker product between the
matrices X and Y.

If TRANS = 'C', y in the conjugate transposed system Z**H *y = scale*b
is solved for, which is equivalent to solve for R and L in

A**H * R + D**H * L = scale * C           (3)
R * B**H + L * E**H = scale * -F

This case (TRANS = 'C') is used to compute an one-norm-based estimate
of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
and (B,E), using ZLACON.

If IJOB >= 1, ZTGSYL computes a Frobenius norm-based estimate of
Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
reciprocal of the smallest singular value of Z.

This is a level-3 BLAS algorithm.
Parameters
 [in] TRANS TRANS is CHARACTER*1 = 'N': solve the generalized sylvester equation (1). = 'C': solve the "conjugate transposed" system (3). [in] IJOB IJOB is INTEGER Specifies what kind of functionality to be performed. =0: solve (1) only. =1: The functionality of 0 and 3. =2: The functionality of 0 and 4. =3: Only an estimate of Dif[(A,D), (B,E)] is computed. (look ahead strategy is used). =4: Only an estimate of Dif[(A,D), (B,E)] is computed. (ZGECON on sub-systems is used). Not referenced if TRANS = 'C'. [in] M M is INTEGER The order of the matrices A and D, and the row dimension of the matrices C, F, R and L. [in] N N is INTEGER The order of the matrices B and E, and the column dimension of the matrices C, F, R and L. [in] A A is COMPLEX*16 array, dimension (LDA, M) The upper triangular matrix A. [in] LDA LDA is INTEGER The leading dimension of the array A. LDA >= max(1, M). [in] B B is COMPLEX*16 array, dimension (LDB, N) The upper triangular matrix B. [in] LDB LDB is INTEGER The leading dimension of the array B. LDB >= max(1, N). [in,out] C C is COMPLEX*16 array, dimension (LDC, N) On entry, C contains the right-hand-side of the first matrix equation in (1) or (3). On exit, if IJOB = 0, 1 or 2, C has been overwritten by the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R, the solution achieved during the computation of the Dif-estimate. [in] LDC LDC is INTEGER The leading dimension of the array C. LDC >= max(1, M). [in] D D is COMPLEX*16 array, dimension (LDD, M) The upper triangular matrix D. [in] LDD LDD is INTEGER The leading dimension of the array D. LDD >= max(1, M). [in] E E is COMPLEX*16 array, dimension (LDE, N) The upper triangular matrix E. [in] LDE LDE is INTEGER The leading dimension of the array E. LDE >= max(1, N). [in,out] F F is COMPLEX*16 array, dimension (LDF, N) On entry, F contains the right-hand-side of the second matrix equation in (1) or (3). On exit, if IJOB = 0, 1 or 2, F has been overwritten by the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L, the solution achieved during the computation of the Dif-estimate. [in] LDF LDF is INTEGER The leading dimension of the array F. LDF >= max(1, M). [out] DIF DIF is DOUBLE PRECISION On exit DIF is the reciprocal of a lower bound of the reciprocal of the Dif-function, i.e. DIF is an upper bound of Dif[(A,D), (B,E)] = sigma-min(Z), where Z as in (2). IF IJOB = 0 or TRANS = 'C', DIF is not referenced. [out] SCALE SCALE is DOUBLE PRECISION On exit SCALE is the scaling factor in (1) or (3). If 0 < SCALE < 1, C and F hold the solutions R and L, resp., to a slightly perturbed system but the input matrices A, B, D and E have not been changed. If SCALE = 0, R and L will hold the solutions to the homogenious system with C = F = 0. [out] WORK WORK is COMPLEX*16 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. LWORK > = 1. If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N). 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] IWORK IWORK is INTEGER array, dimension (M+N+2) [out] INFO INFO is INTEGER =0: successful exit <0: If INFO = -i, the i-th argument had an illegal value. >0: (A, D) and (B, E) have common or very close eigenvalues.
Date
November 2011
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.
References:
 B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software for Solving the Generalized Sylvester Equation and Estimating the Separation between Regular Matrix Pairs, Report UMINF - 93.23, Department of Computing Science, Umea University, S-901 87 Umea, Sweden, December 1993, Revised April 1994, Also as LAPACK Working Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
 B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal. Appl., 15(4):1045-1060, 1994.
 B. Kagstrom and L. Westin, Generalized Schur Methods with Condition Estimators for Solving the Generalized Sylvester Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.

Definition at line 297 of file ztgsyl.f.

297 *
298 * -- LAPACK computational routine (version 3.4.0) --
299 * -- LAPACK is a software package provided by Univ. of Tennessee, --
300 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
301 * November 2011
302 *
303 * .. Scalar Arguments ..
304  CHARACTER trans
305  INTEGER ijob, info, lda, ldb, ldc, ldd, lde, ldf,
306  \$ lwork, m, n
307  DOUBLE PRECISION dif, scale
308 * ..
309 * .. Array Arguments ..
310  INTEGER iwork( * )
311  COMPLEX*16 a( lda, * ), b( ldb, * ), c( ldc, * ),
312  \$ d( ldd, * ), e( lde, * ), f( ldf, * ),
313  \$ work( * )
314 * ..
315 *
316 * =====================================================================
317 * Replaced various illegal calls to CCOPY by calls to CLASET.
318 * Sven Hammarling, 1/5/02.
319 *
320 * .. Parameters ..
321  DOUBLE PRECISION zero, one
322  parameter ( zero = 0.0d+0, one = 1.0d+0 )
323  COMPLEX*16 czero
324  parameter ( czero = (0.0d+0, 0.0d+0) )
325 * ..
326 * .. Local Scalars ..
327  LOGICAL lquery, notran
328  INTEGER i, ie, ifunc, iround, is, isolve, j, je, js, k,
329  \$ linfo, lwmin, mb, nb, p, pq, q
330  DOUBLE PRECISION dscale, dsum, scale2, scaloc
331 * ..
332 * .. External Functions ..
333  LOGICAL lsame
334  INTEGER ilaenv
335  EXTERNAL lsame, ilaenv
336 * ..
337 * .. External Subroutines ..
338  EXTERNAL xerbla, zgemm, zlacpy, zlaset, zscal, ztgsy2
339 * ..
340 * .. Intrinsic Functions ..
341  INTRINSIC dble, dcmplx, max, sqrt
342 * ..
343 * .. Executable Statements ..
344 *
345 * Decode and test input parameters
346 *
347  info = 0
348  notran = lsame( trans, 'N' )
349  lquery = ( lwork.EQ.-1 )
350 *
351  IF( .NOT.notran .AND. .NOT.lsame( trans, 'C' ) ) THEN
352  info = -1
353  ELSE IF( notran ) THEN
354  IF( ( ijob.LT.0 ) .OR. ( ijob.GT.4 ) ) THEN
355  info = -2
356  END IF
357  END IF
358  IF( info.EQ.0 ) THEN
359  IF( m.LE.0 ) THEN
360  info = -3
361  ELSE IF( n.LE.0 ) THEN
362  info = -4
363  ELSE IF( lda.LT.max( 1, m ) ) THEN
364  info = -6
365  ELSE IF( ldb.LT.max( 1, n ) ) THEN
366  info = -8
367  ELSE IF( ldc.LT.max( 1, m ) ) THEN
368  info = -10
369  ELSE IF( ldd.LT.max( 1, m ) ) THEN
370  info = -12
371  ELSE IF( lde.LT.max( 1, n ) ) THEN
372  info = -14
373  ELSE IF( ldf.LT.max( 1, m ) ) THEN
374  info = -16
375  END IF
376  END IF
377 *
378  IF( info.EQ.0 ) THEN
379  IF( notran ) THEN
380  IF( ijob.EQ.1 .OR. ijob.EQ.2 ) THEN
381  lwmin = max( 1, 2*m*n )
382  ELSE
383  lwmin = 1
384  END IF
385  ELSE
386  lwmin = 1
387  END IF
388  work( 1 ) = lwmin
389 *
390  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
391  info = -20
392  END IF
393  END IF
394 *
395  IF( info.NE.0 ) THEN
396  CALL xerbla( 'ZTGSYL', -info )
397  RETURN
398  ELSE IF( lquery ) THEN
399  RETURN
400  END IF
401 *
402 * Quick return if possible
403 *
404  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
405  scale = 1
406  IF( notran ) THEN
407  IF( ijob.NE.0 ) THEN
408  dif = 0
409  END IF
410  END IF
411  RETURN
412  END IF
413 *
414 * Determine optimal block sizes MB and NB
415 *
416  mb = ilaenv( 2, 'ZTGSYL', trans, m, n, -1, -1 )
417  nb = ilaenv( 5, 'ZTGSYL', trans, m, n, -1, -1 )
418 *
419  isolve = 1
420  ifunc = 0
421  IF( notran ) THEN
422  IF( ijob.GE.3 ) THEN
423  ifunc = ijob - 2
424  CALL zlaset( 'F', m, n, czero, czero, c, ldc )
425  CALL zlaset( 'F', m, n, czero, czero, f, ldf )
426  ELSE IF( ijob.GE.1 .AND. notran ) THEN
427  isolve = 2
428  END IF
429  END IF
430 *
431  IF( ( mb.LE.1 .AND. nb.LE.1 ) .OR. ( mb.GE.m .AND. nb.GE.n ) )
432  \$ THEN
433 *
434 * Use unblocked Level 2 solver
435 *
436  DO 30 iround = 1, isolve
437 *
438  scale = one
439  dscale = zero
440  dsum = one
441  pq = m*n
442  CALL ztgsy2( trans, ifunc, m, n, a, lda, b, ldb, c, ldc, d,
443  \$ ldd, e, lde, f, ldf, scale, dsum, dscale,
444  \$ info )
445  IF( dscale.NE.zero ) THEN
446  IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
447  dif = sqrt( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
448  ELSE
449  dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
450  END IF
451  END IF
452  IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
453  IF( notran ) THEN
454  ifunc = ijob
455  END IF
456  scale2 = scale
457  CALL zlacpy( 'F', m, n, c, ldc, work, m )
458  CALL zlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
459  CALL zlaset( 'F', m, n, czero, czero, c, ldc )
460  CALL zlaset( 'F', m, n, czero, czero, f, ldf )
461  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
462  CALL zlacpy( 'F', m, n, work, m, c, ldc )
463  CALL zlacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
464  scale = scale2
465  END IF
466  30 CONTINUE
467 *
468  RETURN
469 *
470  END IF
471 *
472 * Determine block structure of A
473 *
474  p = 0
475  i = 1
476  40 CONTINUE
477  IF( i.GT.m )
478  \$ GO TO 50
479  p = p + 1
480  iwork( p ) = i
481  i = i + mb
482  IF( i.GE.m )
483  \$ GO TO 50
484  GO TO 40
485  50 CONTINUE
486  iwork( p+1 ) = m + 1
487  IF( iwork( p ).EQ.iwork( p+1 ) )
488  \$ p = p - 1
489 *
490 * Determine block structure of B
491 *
492  q = p + 1
493  j = 1
494  60 CONTINUE
495  IF( j.GT.n )
496  \$ GO TO 70
497 *
498  q = q + 1
499  iwork( q ) = j
500  j = j + nb
501  IF( j.GE.n )
502  \$ GO TO 70
503  GO TO 60
504 *
505  70 CONTINUE
506  iwork( q+1 ) = n + 1
507  IF( iwork( q ).EQ.iwork( q+1 ) )
508  \$ q = q - 1
509 *
510  IF( notran ) THEN
511  DO 150 iround = 1, isolve
512 *
513 * Solve (I, J) - subsystem
514 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
515 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
516 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
517 *
518  pq = 0
519  scale = one
520  dscale = zero
521  dsum = one
522  DO 130 j = p + 2, q
523  js = iwork( j )
524  je = iwork( j+1 ) - 1
525  nb = je - js + 1
526  DO 120 i = p, 1, -1
527  is = iwork( i )
528  ie = iwork( i+1 ) - 1
529  mb = ie - is + 1
530  CALL ztgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
531  \$ b( js, js ), ldb, c( is, js ), ldc,
532  \$ d( is, is ), ldd, e( js, js ), lde,
533  \$ f( is, js ), ldf, scaloc, dsum, dscale,
534  \$ linfo )
535  IF( linfo.GT.0 )
536  \$ info = linfo
537  pq = pq + mb*nb
538  IF( scaloc.NE.one ) THEN
539  DO 80 k = 1, js - 1
540  CALL zscal( m, dcmplx( scaloc, zero ),
541  \$ c( 1, k ), 1 )
542  CALL zscal( m, dcmplx( scaloc, zero ),
543  \$ f( 1, k ), 1 )
544  80 CONTINUE
545  DO 90 k = js, je
546  CALL zscal( is-1, dcmplx( scaloc, zero ),
547  \$ c( 1, k ), 1 )
548  CALL zscal( is-1, dcmplx( scaloc, zero ),
549  \$ f( 1, k ), 1 )
550  90 CONTINUE
551  DO 100 k = js, je
552  CALL zscal( m-ie, dcmplx( scaloc, zero ),
553  \$ c( ie+1, k ), 1 )
554  CALL zscal( m-ie, dcmplx( scaloc, zero ),
555  \$ f( ie+1, k ), 1 )
556  100 CONTINUE
557  DO 110 k = je + 1, n
558  CALL zscal( m, dcmplx( scaloc, zero ),
559  \$ c( 1, k ), 1 )
560  CALL zscal( m, dcmplx( scaloc, zero ),
561  \$ f( 1, k ), 1 )
562  110 CONTINUE
563  scale = scale*scaloc
564  END IF
565 *
566 * Substitute R(I,J) and L(I,J) into remaining equation.
567 *
568  IF( i.GT.1 ) THEN
569  CALL zgemm( 'N', 'N', is-1, nb, mb,
570  \$ dcmplx( -one, zero ), a( 1, is ), lda,
571  \$ c( is, js ), ldc, dcmplx( one, zero ),
572  \$ c( 1, js ), ldc )
573  CALL zgemm( 'N', 'N', is-1, nb, mb,
574  \$ dcmplx( -one, zero ), d( 1, is ), ldd,
575  \$ c( is, js ), ldc, dcmplx( one, zero ),
576  \$ f( 1, js ), ldf )
577  END IF
578  IF( j.LT.q ) THEN
579  CALL zgemm( 'N', 'N', mb, n-je, nb,
580  \$ dcmplx( one, zero ), f( is, js ), ldf,
581  \$ b( js, je+1 ), ldb,
582  \$ dcmplx( one, zero ), c( is, je+1 ),
583  \$ ldc )
584  CALL zgemm( 'N', 'N', mb, n-je, nb,
585  \$ dcmplx( one, zero ), f( is, js ), ldf,
586  \$ e( js, je+1 ), lde,
587  \$ dcmplx( one, zero ), f( is, je+1 ),
588  \$ ldf )
589  END IF
590  120 CONTINUE
591  130 CONTINUE
592  IF( dscale.NE.zero ) THEN
593  IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
594  dif = sqrt( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
595  ELSE
596  dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
597  END IF
598  END IF
599  IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
600  IF( notran ) THEN
601  ifunc = ijob
602  END IF
603  scale2 = scale
604  CALL zlacpy( 'F', m, n, c, ldc, work, m )
605  CALL zlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
606  CALL zlaset( 'F', m, n, czero, czero, c, ldc )
607  CALL zlaset( 'F', m, n, czero, czero, f, ldf )
608  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
609  CALL zlacpy( 'F', m, n, work, m, c, ldc )
610  CALL zlacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
611  scale = scale2
612  END IF
613  150 CONTINUE
614  ELSE
615 *
616 * Solve transposed (I, J)-subsystem
617 * A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J)
618 * R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
619 * for I = 1,2,..., P; J = Q, Q-1,..., 1
620 *
621  scale = one
622  DO 210 i = 1, p
623  is = iwork( i )
624  ie = iwork( i+1 ) - 1
625  mb = ie - is + 1
626  DO 200 j = q, p + 2, -1
627  js = iwork( j )
628  je = iwork( j+1 ) - 1
629  nb = je - js + 1
630  CALL ztgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
631  \$ b( js, js ), ldb, c( is, js ), ldc,
632  \$ d( is, is ), ldd, e( js, js ), lde,
633  \$ f( is, js ), ldf, scaloc, dsum, dscale,
634  \$ linfo )
635  IF( linfo.GT.0 )
636  \$ info = linfo
637  IF( scaloc.NE.one ) THEN
638  DO 160 k = 1, js - 1
639  CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
640  \$ 1 )
641  CALL zscal( m, dcmplx( scaloc, zero ), f( 1, k ),
642  \$ 1 )
643  160 CONTINUE
644  DO 170 k = js, je
645  CALL zscal( is-1, dcmplx( scaloc, zero ),
646  \$ c( 1, k ), 1 )
647  CALL zscal( is-1, dcmplx( scaloc, zero ),
648  \$ f( 1, k ), 1 )
649  170 CONTINUE
650  DO 180 k = js, je
651  CALL zscal( m-ie, dcmplx( scaloc, zero ),
652  \$ c( ie+1, k ), 1 )
653  CALL zscal( m-ie, dcmplx( scaloc, zero ),
654  \$ f( ie+1, k ), 1 )
655  180 CONTINUE
656  DO 190 k = je + 1, n
657  CALL zscal( m, dcmplx( scaloc, zero ), c( 1, k ),
658  \$ 1 )
659  CALL zscal( m, dcmplx( scaloc, zero ), f( 1, k ),
660  \$ 1 )
661  190 CONTINUE
662  scale = scale*scaloc
663  END IF
664 *
665 * Substitute R(I,J) and L(I,J) into remaining equation.
666 *
667  IF( j.GT.p+2 ) THEN
668  CALL zgemm( 'N', 'C', mb, js-1, nb,
669  \$ dcmplx( one, zero ), c( is, js ), ldc,
670  \$ b( 1, js ), ldb, dcmplx( one, zero ),
671  \$ f( is, 1 ), ldf )
672  CALL zgemm( 'N', 'C', mb, js-1, nb,
673  \$ dcmplx( one, zero ), f( is, js ), ldf,
674  \$ e( 1, js ), lde, dcmplx( one, zero ),
675  \$ f( is, 1 ), ldf )
676  END IF
677  IF( i.LT.p ) THEN
678  CALL zgemm( 'C', 'N', m-ie, nb, mb,
679  \$ dcmplx( -one, zero ), a( is, ie+1 ), lda,
680  \$ c( is, js ), ldc, dcmplx( one, zero ),
681  \$ c( ie+1, js ), ldc )
682  CALL zgemm( 'C', 'N', m-ie, nb, mb,
683  \$ dcmplx( -one, zero ), d( is, ie+1 ), ldd,
684  \$ f( is, js ), ldf, dcmplx( one, zero ),
685  \$ c( ie+1, js ), ldc )
686  END IF
687  200 CONTINUE
688  210 CONTINUE
689  END IF
690 *
691  work( 1 ) = lwmin
692 *
693  RETURN
694 *
695 * End of ZTGSYL
696 *
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine ztgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, INFO)
ZTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition: ztgsy2.f:261
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: