Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[8425] trunk/octave-forge/main/control/devel
|
From: <par...@us...> - 2011-07-31 19:51:25
|
Revision: 8425
http://octave.svn.sourceforge.net/octave/?rev=8425&view=rev
Author: paramaniac
Date: 2011-07-31 19:51:16 +0000 (Sun, 31 Jul 2011)
Log Message:
-----------
control: add ncfsyn (continuous-time case only)
Added Paths:
-----------
trunk/octave-forge/main/control/devel/MDSSystem.m
trunk/octave-forge/main/control/devel/ncfsyn/
trunk/octave-forge/main/control/devel/ncfsyn/MA02AD.f
trunk/octave-forge/main/control/devel/ncfsyn/MA02ED.f
trunk/octave-forge/main/control/devel/ncfsyn/MA02GD.f
trunk/octave-forge/main/control/devel/ncfsyn/MB01RU.f
trunk/octave-forge/main/control/devel/ncfsyn/MB01RX.f
trunk/octave-forge/main/control/devel/ncfsyn/MB01RY.f
trunk/octave-forge/main/control/devel/ncfsyn/MB01SD.f
trunk/octave-forge/main/control/devel/ncfsyn/MB01UD.f
trunk/octave-forge/main/control/devel/ncfsyn/MB02PD.f
trunk/octave-forge/main/control/devel/ncfsyn/MB02VD.f
trunk/octave-forge/main/control/devel/ncfsyn/SB02MR.f
trunk/octave-forge/main/control/devel/ncfsyn/SB02MS.f
trunk/octave-forge/main/control/devel/ncfsyn/SB02MV.f
trunk/octave-forge/main/control/devel/ncfsyn/SB02MW.f
trunk/octave-forge/main/control/devel/ncfsyn/SB02QD.f
trunk/octave-forge/main/control/devel/ncfsyn/SB02RD.f
trunk/octave-forge/main/control/devel/ncfsyn/SB02RU.f
trunk/octave-forge/main/control/devel/ncfsyn/SB02SD.f
trunk/octave-forge/main/control/devel/ncfsyn/SB03MV.f
trunk/octave-forge/main/control/devel/ncfsyn/SB03MW.f
trunk/octave-forge/main/control/devel/ncfsyn/SB03MX.f
trunk/octave-forge/main/control/devel/ncfsyn/SB03MY.f
trunk/octave-forge/main/control/devel/ncfsyn/SB03QX.f
trunk/octave-forge/main/control/devel/ncfsyn/SB03QY.f
trunk/octave-forge/main/control/devel/ncfsyn/SB03SX.f
trunk/octave-forge/main/control/devel/ncfsyn/SB03SY.f
trunk/octave-forge/main/control/devel/ncfsyn/SB04PX.f
trunk/octave-forge/main/control/devel/ncfsyn/SB10ID.f
trunk/octave-forge/main/control/devel/ncfsyn/SB10JD.f
trunk/octave-forge/main/control/devel/ncfsyn/common.cc
trunk/octave-forge/main/control/devel/ncfsyn/makefile_ncfsyn.m
trunk/octave-forge/main/control/devel/ncfsyn/ncfsyn.m
trunk/octave-forge/main/control/devel/ncfsyn/select.f
trunk/octave-forge/main/control/devel/ncfsyn/slsb10id.cc
Added: trunk/octave-forge/main/control/devel/MDSSystem.m
===================================================================
--- trunk/octave-forge/main/control/devel/MDSSystem.m (rev 0)
+++ trunk/octave-forge/main/control/devel/MDSSystem.m 2011-07-31 19:51:16 UTC (rev 8425)
@@ -0,0 +1,123 @@
+% Tabula Rasa
+clear all, close all, clc
+
+% Nominal Values
+m_nom = 3; % Mass
+c_nom = 1; % Damping Coefficient
+k_nom = 2; % Spring Stiffness
+
+% Perturbations
+p_m = 0.4; % 40% uncertainty in the mass
+p_c = 0.2; % 20% uncertainty in the damping coefficient
+p_k = 0.3; % 30% uncertainty in the spring stiffness
+
+% +---------------+
+% | d_m 0 0 |
+% +-----| 0 d_c 0 |<----+
+% u_m | | 0 0 d_k | | y_m
+% u_c | +---------------+ | y_c
+% u_k | | y_k
+% | +---------------+ |
+% +---->| |-----+
+% | G_nom |
+% u ----->| |-----> y
+% +---------------+
+
+% State-Space Representation
+A = [ 0, 1
+ -k_nom/m_nom, -c_nom/m_nom ];
+
+B1 = [ 0, 0, 0
+ -p_m, -p_c/m_nom, -p_k/m_nom ];
+
+B2 = [ 0
+ 1/m_nom ];
+
+C1 = [ -k_nom/m_nom, -c_nom/m_nom
+ 0, c_nom
+ k_nom, 0 ];
+
+C2 = [ 1, 0 ];
+
+D11 = [ -p_m, -p_c/m_nom, -p_k/m_nom
+ 0, 0, 0
+ 0, 0, 0 ];
+
+D12 = [ 1/m_nom
+ 0
+ 0 ];
+
+D21 = [ 0, 0, 0 ];
+
+D22 = [ 0 ];
+
+inname = {"u_m", "u_c", "u_k", "u"}; % Input Names
+outname = {"y_m", "y_c", "y_k", "y"}; % Output Names
+
+G_nom = ss (A, [B1, B2], [C1; C2], [D11, D12; D21, D22], \
+ "inname", inname, "outname", outname);
+
+% +-------+
+% +--------------------->| W_p |----------> e_p
+% | +-------+
+% | +-------+
+% | +---->| W_u |----------> e_u
+% | | +-------+
+% | | +---------+
+% | | ->| |->
+% r + e | +-------+ u | | G_nom |
+% ----->(+)---+-->| K |----+--->| |----+----> y
+% ^ - +-------+ +---------+ |
+% | |
+% +-----------------------------------------+
+
+% Weighting Functions
+s = tf ("s");
+W_p = 0.95 * (s^2 + 1.8*s + 10) / (s^2 + 8.0*s + 0.01);
+W_u = 10^-2;
+
+% Suboptimal H-infinity Controller Design
+G = G_nom(4, 4); % Extract output y and input u
+%K = mixsyn (G, W_p, W_u);
+%[K, ~, gamma] = mixsyn (G, W_p, W_u)
+[K, N, gamma] = mixsyn (G, W_p, W_u)
+% g_max = 10;
+% K = mixsyn (G, W_p, W_u, [], g_max)
+% [K, ~, gamma] = mixsyn (G, W_p, W_u, [], g_max)
+
+% Open Loop
+L = G * K;
+
+% Closed Loop
+T = feedback (L);
+
+figure (1)
+sigma (T)
+
+figure (2)
+sigma (N)
+
+norm (N, inf)
+
+figure (3)
+step (T)
+
+
+figure (4)
+
+w = logspace (-1, 1, 100);
+
+% Uncertainties
+% -1 <= delta_m, delta_c, delta_k <= 1
+[delta_m, delta_c, delta_k] = ndgrid ([-1, 0, 1], [-1, 0, 1], [-1, 0, 1]);
+
+for k = 1 : numel (delta_1)
+ Delta = diag ([delta_m(k), delta_c(k), delta_k(k)]);
+ G_per = lft (Delta, G_nom);
+ %figure (4)
+ bode (G_per, w)
+ subplot (2, 1, 1)
+ hold on
+ subplot (2, 1, 2)
+ hold on
+endfor
\ No newline at end of file
Property changes on: trunk/octave-forge/main/control/devel/ncfsyn
___________________________________________________________________
Added: svn:ignore
+ *.oct
Added: trunk/octave-forge/main/control/devel/ncfsyn/MA02AD.f
===================================================================
--- trunk/octave-forge/main/control/devel/ncfsyn/MA02AD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ncfsyn/MA02AD.f 2011-07-31 19:51:16 UTC (rev 8425)
@@ -0,0 +1,108 @@
+ SUBROUTINE MA02AD( JOB, M, N, A, LDA, B, LDB )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2010 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To transpose all or part of a two-dimensional matrix A into
+C another matrix B.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C JOB CHARACTER*1
+C Specifies the part of the matrix A to be transposed into B
+C as follows:
+C = 'U': Upper triangular part;
+C = 'L': Lower triangular part;
+C Otherwise: All of the matrix A.
+C
+C Input/Output Parameters
+C
+C M (input) INTEGER
+C The number of rows of the matrix A. M >= 0.
+C
+C N (input) INTEGER
+C The number of columns of the matrix A. N >= 0.
+C
+C A (input) DOUBLE PRECISION array, dimension (LDA,N)
+C The m-by-n matrix A. If JOB = 'U', only the upper
+C triangle or trapezoid is accessed; if JOB = 'L', only the
+C lower triangle or trapezoid is accessed.
+C
+C LDA INTEGER
+C The leading dimension of the array A. LDA >= max(1,M).
+C
+C B (output) DOUBLE PRECISION array, dimension (LDB,M)
+C B = A' in the locations specified by JOB.
+C
+C LDB INTEGER
+C The leading dimension of the array B. LDB >= max(1,N).
+C
+C CONTRIBUTOR
+C
+C A. Varga, German Aerospace Center,
+C DLR Oberpfaffenhofen, March 1998.
+C Based on the RASP routine DMTRA.
+C
+C REVISIONS
+C
+C -
+C
+C ******************************************************************
+C
+C .. Scalar Arguments ..
+ CHARACTER JOB
+ INTEGER LDA, LDB, M, N
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*), B(LDB,*)
+C .. Local Scalars ..
+ INTEGER I, J
+C .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C .. Intrinsic Functions ..
+ INTRINSIC MIN
+C
+C .. Executable Statements ..
+C
+ IF( LSAME( JOB, 'U' ) ) THEN
+ DO 20 J = 1, N
+ DO 10 I = 1, MIN( J, M )
+ B(J,I) = A(I,J)
+ 10 CONTINUE
+ 20 CONTINUE
+ ELSE IF( LSAME( JOB, 'L' ) ) THEN
+ DO 40 J = 1, N
+ DO 30 I = J, M
+ B(J,I) = A(I,J)
+ 30 CONTINUE
+ 40 CONTINUE
+ ELSE
+ DO 60 J = 1, N
+ DO 50 I = 1, M
+ B(J,I) = A(I,J)
+ 50 CONTINUE
+ 60 CONTINUE
+ END IF
+C
+ RETURN
+C *** Last line of MA02AD ***
+ END
Property changes on: trunk/octave-forge/main/control/devel/ncfsyn/MA02AD.f
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/control/devel/ncfsyn/MA02ED.f
===================================================================
--- trunk/octave-forge/main/control/devel/ncfsyn/MA02ED.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ncfsyn/MA02ED.f 2011-07-31 19:51:16 UTC (rev 8425)
@@ -0,0 +1,99 @@
+ SUBROUTINE MA02ED( UPLO, N, A, LDA )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2010 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To store by symmetry the upper or lower triangle of a symmetric
+C matrix, given the other triangle.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C UPLO CHARACTER*1
+C Specifies which part of the matrix is given as follows:
+C = 'U': Upper triangular part;
+C = 'L': Lower triangular part.
+C For all other values, the array A is not referenced.
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The order of the matrix A. N >= 0.
+C
+C A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+C On entry, the leading N-by-N upper triangular part
+C (if UPLO = 'U'), or lower triangular part (if UPLO = 'L'),
+C of this array must contain the corresponding upper or
+C lower triangle of the symmetric matrix A.
+C On exit, the leading N-by-N part of this array contains
+C the symmetric matrix A with all elements stored.
+C
+C LDA INTEGER
+C The leading dimension of the array A. LDA >= max(1,N).
+C
+C CONTRIBUTOR
+C
+C V. Sima, Research Institute for Informatics, Bucharest, Romania,
+C Oct. 1998.
+C
+C REVISIONS
+C
+C -
+C
+C ******************************************************************
+C
+C .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER LDA, N
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*)
+C .. Local Scalars ..
+ INTEGER J
+C .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C .. External Subroutines ..
+ EXTERNAL DCOPY
+C
+C .. Executable Statements ..
+C
+C For efficiency reasons, the parameters are not checked for errors.
+C
+ IF( LSAME( UPLO, 'L' ) ) THEN
+C
+C Construct the upper triangle of A.
+C
+ DO 20 J = 2, N
+ CALL DCOPY( J-1, A(J,1), LDA, A(1,J), 1 )
+ 20 CONTINUE
+C
+ ELSE IF( LSAME( UPLO, 'U' ) ) THEN
+C
+C Construct the lower triangle of A.
+C
+ DO 40 J = 2, N
+ CALL DCOPY( J-1, A(1,J), 1, A(J,1), LDA )
+ 40 CONTINUE
+C
+ END IF
+ RETURN
+C *** Last line of MA02ED ***
+ END
Property changes on: trunk/octave-forge/main/control/devel/ncfsyn/MA02ED.f
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/control/devel/ncfsyn/MA02GD.f
===================================================================
--- trunk/octave-forge/main/control/devel/ncfsyn/MA02GD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ncfsyn/MA02GD.f 2011-07-31 19:51:16 UTC (rev 8425)
@@ -0,0 +1,158 @@
+ SUBROUTINE MA02GD( N, A, LDA, K1, K2, IPIV, INCX )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2010 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To perform a series of column interchanges on the matrix A.
+C One column interchange is initiated for each of columns K1 through
+C K2 of A. This is useful for solving linear systems X*A = B, when
+C the matrix A has already been factored by LAPACK Library routine
+C DGETRF.
+C
+C ARGUMENTS
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The number of rows of the matrix A. N >= 0.
+C
+C A (input/output) DOUBLE PRECISION array, dimension (LDA,*)
+C On entry, the leading N-by-M part of this array must
+C contain the matrix A to which the column interchanges will
+C be applied, where M is the largest element of IPIV(K), for
+C K = K1, ..., K2.
+C On exit, the leading N-by-M part of this array contains
+C the permuted matrix.
+C
+C LDA INTEGER
+C The leading dimension of the array A. LDA >= MAX(1,N).
+C
+C K1 (input) INTEGER
+C The first element of IPIV for which a column interchange
+C will be done.
+C
+C K2 (input) INTEGER
+C The last element of IPIV for which a column interchange
+C will be done.
+C
+C IPIV (input) INTEGER array, dimension (K1+(K2-K1)*abs(INCX))
+C The vector of interchanging (pivot) indices. Only the
+C elements in positions K1 through K2 of IPIV are accessed.
+C IPIV(K) = L implies columns K and L are to be
+C interchanged.
+C
+C INCX (input) INTEGER
+C The increment between successive values of IPIV.
+C If INCX is negative, the interchanges are applied in
+C reverse order.
+C
+C METHOD
+C
+C The columns IPIV(K) and K are swapped for K = K1, ..., K2, for
+C INCX = 1 (and similarly, for INCX <> 1).
+C
+C FURTHER COMMENTS
+C
+C This routine is the column-oriented counterpart of the LAPACK
+C Library routine DLASWP. The LAPACK Library routine DLAPMT cannot
+C be used in this context. To solve the system X*A = B, where A and
+C B are N-by-N and M-by-N, respectively, the following statements
+C can be used:
+C
+C CALL DGETRF( N, N, A, LDA, IPIV, INFO )
+C CALL DTRSM( 'R', 'U', 'N', 'N', M, N, ONE, A, LDA, B, LDB )
+C CALL DTRSM( 'R', 'L', 'N', 'U', M, N, ONE, A, LDA, B, LDB )
+C CALL MA02GD( M, B, LDB, 1, N, IPIV, -1 )
+C
+C CONTRIBUTOR
+C
+C V. Sima, Research Institute for Informatics, Bucharest, Mar. 2000.
+C
+C REVISIONS
+C
+C V. Sima, Research Institute for Informatics, Bucharest, Apr. 2008.
+C
+C KEYWORDS
+C
+C Elementary matrix operations, linear algebra.
+C
+C ******************************************************************
+C
+C .. Scalar Arguments ..
+ INTEGER INCX, K1, K2, LDA, N
+C ..
+C .. Array Arguments ..
+ INTEGER IPIV( * )
+ DOUBLE PRECISION A( LDA, * )
+C ..
+C .. Local Scalars ..
+ INTEGER J, JP, JX
+C ..
+C .. External Subroutines ..
+ EXTERNAL DSWAP
+C ..
+C .. Executable Statements ..
+C
+C Quick return if possible.
+C
+ IF( INCX.EQ.0 .OR. N.EQ.0 )
+ $ RETURN
+C
+C Interchange column J with column IPIV(J) for each of columns K1
+C through K2.
+C
+ IF( INCX.GT.0 ) THEN
+ JX = K1
+ ELSE
+ JX = 1 + ( 1-K2 )*INCX
+ END IF
+C
+ IF( INCX.EQ.1 ) THEN
+C
+ DO 10 J = K1, K2
+ JP = IPIV( J )
+ IF( JP.NE.J )
+ $ CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
+ 10 CONTINUE
+C
+ ELSE IF( INCX.GT.1 ) THEN
+C
+ DO 20 J = K1, K2
+ JP = IPIV( JX )
+ IF( JP.NE.J )
+ $ CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
+ JX = JX + INCX
+ 20 CONTINUE
+C
+ ELSE IF( INCX.LT.0 ) THEN
+C
+ DO 30 J = K2, K1, -1
+ JP = IPIV( JX )
+ IF( JP.NE.J )
+ $ CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
+ JX = JX + INCX
+ 30 CONTINUE
+C
+ END IF
+C
+ RETURN
+C
+C *** Last line of MA02GD ***
+ END
Property changes on: trunk/octave-forge/main/control/devel/ncfsyn/MA02GD.f
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/control/devel/ncfsyn/MB01RU.f
===================================================================
--- trunk/octave-forge/main/control/devel/ncfsyn/MB01RU.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ncfsyn/MB01RU.f 2011-07-31 19:51:16 UTC (rev 8425)
@@ -0,0 +1,282 @@
+ SUBROUTINE MB01RU( UPLO, TRANS, M, N, ALPHA, BETA, R, LDR, A, LDA,
+ $ X, LDX, DWORK, LDWORK, INFO )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2010 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To compute the matrix formula
+C _
+C R = alpha*R + beta*op( A )*X*op( A )',
+C _
+C where alpha and beta are scalars, R, X, and R are symmetric
+C matrices, A is a general matrix, and op( A ) is one of
+C
+C op( A ) = A or op( A ) = A'.
+C
+C The result is overwritten on R.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C UPLO CHARACTER*1
+C Specifies which triangles of the symmetric matrices R
+C and X are given as follows:
+C = 'U': the upper triangular part is given;
+C = 'L': the lower triangular part is given.
+C
+C TRANS CHARACTER*1
+C Specifies the form of op( A ) to be used in the matrix
+C multiplication as follows:
+C = 'N': op( A ) = A;
+C = 'T': op( A ) = A';
+C = 'C': op( A ) = A'.
+C
+C Input/Output Parameters
+C
+C M (input) INTEGER _
+C The order of the matrices R and R and the number of rows
+C of the matrix op( A ). M >= 0.
+C
+C N (input) INTEGER
+C The order of the matrix X and the number of columns of the
+C the matrix op( A ). N >= 0.
+C
+C ALPHA (input) DOUBLE PRECISION
+C The scalar alpha. When alpha is zero then R need not be
+C set before entry, except when R is identified with X in
+C the call.
+C
+C BETA (input) DOUBLE PRECISION
+C The scalar beta. When beta is zero then A and X are not
+C referenced.
+C
+C R (input/output) DOUBLE PRECISION array, dimension (LDR,M)
+C On entry with UPLO = 'U', the leading M-by-M upper
+C triangular part of this array must contain the upper
+C triangular part of the symmetric matrix R.
+C On entry with UPLO = 'L', the leading M-by-M lower
+C triangular part of this array must contain the lower
+C triangular part of the symmetric matrix R.
+C On exit, the leading M-by-M upper triangular part (if
+C UPLO = 'U'), or lower triangular part (if UPLO = 'L'), of
+C this array contains the corresponding triangular part of
+C _
+C the computed matrix R.
+C
+C LDR INTEGER
+C The leading dimension of array R. LDR >= MAX(1,M).
+C
+C A (input) DOUBLE PRECISION array, dimension (LDA,k)
+C where k is N when TRANS = 'N' and is M when TRANS = 'T' or
+C TRANS = 'C'.
+C On entry with TRANS = 'N', the leading M-by-N part of this
+C array must contain the matrix A.
+C On entry with TRANS = 'T' or TRANS = 'C', the leading
+C N-by-M part of this array must contain the matrix A.
+C
+C LDA INTEGER
+C The leading dimension of array A. LDA >= MAX(1,k),
+C where k is M when TRANS = 'N' and is N when TRANS = 'T' or
+C TRANS = 'C'.
+C
+C X (input) DOUBLE PRECISION array, dimension (LDX,N)
+C On entry, if UPLO = 'U', the leading N-by-N upper
+C triangular part of this array must contain the upper
+C triangular part of the symmetric matrix X and the strictly
+C lower triangular part of the array is not referenced.
+C On entry, if UPLO = 'L', the leading N-by-N lower
+C triangular part of this array must contain the lower
+C triangular part of the symmetric matrix X and the strictly
+C upper triangular part of the array is not referenced.
+C The diagonal elements of this array are modified
+C internally, but are restored on exit.
+C
+C LDX INTEGER
+C The leading dimension of array X. LDX >= MAX(1,N).
+C
+C Workspace
+C
+C DWORK DOUBLE PRECISION array, dimension (LDWORK)
+C This array is not referenced when beta = 0, or M*N = 0.
+C
+C LDWORK The length of the array DWORK.
+C LDWORK >= M*N, if beta <> 0;
+C LDWORK >= 0, if beta = 0.
+C
+C Error Indicator
+C
+C INFO INTEGER
+C = 0: successful exit;
+C < 0: if INFO = -k, the k-th argument had an illegal
+C value.
+C
+C METHOD
+C
+C The matrix expression is efficiently evaluated taking the symmetry
+C into account. Specifically, let X = T + T', with T an upper or
+C lower triangular matrix, defined by
+C
+C T = triu( X ) - (1/2)*diag( X ), if UPLO = 'U',
+C T = tril( X ) - (1/2)*diag( X ), if UPLO = 'L',
+C
+C where triu, tril, and diag denote the upper triangular part, lower
+C triangular part, and diagonal part of X, respectively. Then,
+C
+C A*X*A' = ( A*T )*A' + A*( A*T )', for TRANS = 'N',
+C A'*X*A = A'*( T*A ) + ( T*A )'*A, for TRANS = 'T', or 'C',
+C
+C which involve BLAS 3 operations (DTRMM and DSYR2K).
+C
+C NUMERICAL ASPECTS
+C
+C The algorithm requires approximately
+C
+C 2 2
+C 3/2 x M x N + 1/2 x M
+C
+C operations.
+C
+C FURTHER COMMENTS
+C
+C This is a simpler version for MB01RD.
+C
+C CONTRIBUTORS
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, Jan. 1999.
+C
+C REVISIONS
+C
+C A. Varga, German Aerospace Center, Oberpfaffenhofen, March 2004.
+C V. Sima, Research Institute for Informatics, Bucharest, Mar. 2004.
+C
+C KEYWORDS
+C
+C Elementary matrix operations, matrix algebra, matrix operations.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE, TWO, HALF
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
+ $ HALF = 0.5D0 )
+C .. Scalar Arguments ..
+ CHARACTER TRANS, UPLO
+ INTEGER INFO, LDA, LDR, LDWORK, LDX, M, N
+ DOUBLE PRECISION ALPHA, BETA
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*), DWORK(*), R(LDR,*), X(LDX,*)
+C .. Local Scalars ..
+ LOGICAL LTRANS, LUPLO
+C .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C .. External Subroutines ..
+ EXTERNAL DLACPY, DLASCL, DLASET, DSCAL, DSYR2K, DTRMM,
+ $ XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC MAX
+C .. Executable Statements ..
+C
+C Test the input scalar arguments.
+C
+ INFO = 0
+ LUPLO = LSAME( UPLO, 'U' )
+ LTRANS = LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' )
+C
+ IF( ( .NOT.LUPLO ).AND.( .NOT.LSAME( UPLO, 'L' ) ) )THEN
+ INFO = -1
+ ELSE IF( ( .NOT.LTRANS ).AND.( .NOT.LSAME( TRANS, 'N' ) ) )THEN
+ INFO = -2
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDR.LT.MAX( 1, M ) ) THEN
+ INFO = -8
+ ELSE IF( LDA.LT.1 .OR. ( LTRANS .AND. LDA.LT.N ) .OR.
+ $ ( .NOT.LTRANS .AND. LDA.LT.M ) ) THEN
+ INFO = -10
+ ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
+ INFO = -12
+ ELSE IF( ( BETA.NE.ZERO .AND. LDWORK.LT.M*N )
+ $ .OR.( BETA.EQ.ZERO .AND. LDWORK.LT.0 ) ) THEN
+ INFO = -14
+ END IF
+C
+ IF ( INFO.NE.0 ) THEN
+C
+C Error return.
+C
+ CALL XERBLA( 'MB01RU', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ IF ( M.EQ.0 )
+ $ RETURN
+C
+ IF ( BETA.EQ.ZERO .OR. N.EQ.0 ) THEN
+ IF ( ALPHA.EQ.ZERO ) THEN
+C
+C Special case alpha = 0.
+C
+ CALL DLASET( UPLO, M, M, ZERO, ZERO, R, LDR )
+ ELSE
+C
+C Special case beta = 0 or N = 0.
+C
+ IF ( ALPHA.NE.ONE )
+ $ CALL DLASCL( UPLO, 0, 0, ONE, ALPHA, M, M, R, LDR, INFO )
+ END IF
+ RETURN
+ END IF
+C
+C General case: beta <> 0.
+C Compute W = op( A )*T or W = T*op( A ) in DWORK, and apply the
+C updating formula (see METHOD section).
+C Workspace: need M*N.
+C
+ CALL DSCAL( N, HALF, X, LDX+1 )
+C
+ IF( LTRANS ) THEN
+C
+ CALL DLACPY( 'Full', N, M, A, LDA, DWORK, N )
+ CALL DTRMM( 'Left', UPLO, 'NoTranspose', 'Non-unit', N, M,
+ $ ONE, X, LDX, DWORK, N )
+ CALL DSYR2K( UPLO, TRANS, M, N, BETA, DWORK, N, A, LDA, ALPHA,
+ $ R, LDR )
+C
+ ELSE
+C
+ CALL DLACPY( 'Full', M, N, A, LDA, DWORK, M )
+ CALL DTRMM( 'Right', UPLO, 'NoTranspose', 'Non-unit', M, N,
+ $ ONE, X, LDX, DWORK, M )
+ CALL DSYR2K( UPLO, TRANS, M, N, BETA, DWORK, M, A, LDA, ALPHA,
+ $ R, LDR )
+C
+ END IF
+C
+ CALL DSCAL( N, TWO, X, LDX+1 )
+C
+ RETURN
+C *** Last line of MB01RU ***
+ END
Property changes on: trunk/octave-forge/main/control/devel/ncfsyn/MB01RU.f
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/control/devel/ncfsyn/MB01RX.f
===================================================================
--- trunk/octave-forge/main/control/devel/ncfsyn/MB01RX.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ncfsyn/MB01RX.f 2011-07-31 19:51:16 UTC (rev 8425)
@@ -0,0 +1,315 @@
+ SUBROUTINE MB01RX( SIDE, UPLO, TRANS, M, N, ALPHA, BETA, R, LDR,
+ $ A, LDA, B, LDB, INFO )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2010 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To compute either the upper or lower triangular part of one of the
+C matrix formulas
+C _
+C R = alpha*R + beta*op( A )*B, (1)
+C _
+C R = alpha*R + beta*B*op( A ), (2)
+C _
+C where alpha and beta are scalars, R and R are m-by-m matrices,
+C op( A ) and B are m-by-n and n-by-m matrices for (1), or n-by-m
+C and m-by-n matrices for (2), respectively, and op( A ) is one of
+C
+C op( A ) = A or op( A ) = A', the transpose of A.
+C
+C The result is overwritten on R.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C SIDE CHARACTER*1
+C Specifies whether the matrix A appears on the left or
+C right in the matrix product as follows:
+C _
+C = 'L': R = alpha*R + beta*op( A )*B;
+C _
+C = 'R': R = alpha*R + beta*B*op( A ).
+C
+C UPLO CHARACTER*1 _
+C Specifies which triangles of the matrices R and R are
+C computed and given, respectively, as follows:
+C = 'U': the upper triangular part;
+C = 'L': the lower triangular part.
+C
+C TRANS CHARACTER*1
+C Specifies the form of op( A ) to be used in the matrix
+C multiplication as follows:
+C = 'N': op( A ) = A;
+C = 'T': op( A ) = A';
+C = 'C': op( A ) = A'.
+C
+C Input/Output Parameters
+C
+C M (input) INTEGER _
+C The order of the matrices R and R, the number of rows of
+C the matrix op( A ) and the number of columns of the
+C matrix B, for SIDE = 'L', or the number of rows of the
+C matrix B and the number of columns of the matrix op( A ),
+C for SIDE = 'R'. M >= 0.
+C
+C N (input) INTEGER
+C The number of rows of the matrix B and the number of
+C columns of the matrix op( A ), for SIDE = 'L', or the
+C number of rows of the matrix op( A ) and the number of
+C columns of the matrix B, for SIDE = 'R'. N >= 0.
+C
+C ALPHA (input) DOUBLE PRECISION
+C The scalar alpha. When alpha is zero then R need not be
+C set before entry.
+C
+C BETA (input) DOUBLE PRECISION
+C The scalar beta. When beta is zero then A and B are not
+C referenced.
+C
+C R (input/output) DOUBLE PRECISION array, dimension (LDR,M)
+C On entry with UPLO = 'U', the leading M-by-M upper
+C triangular part of this array must contain the upper
+C triangular part of the matrix R; the strictly lower
+C triangular part of the array is not referenced.
+C On entry with UPLO = 'L', the leading M-by-M lower
+C triangular part of this array must contain the lower
+C triangular part of the matrix R; the strictly upper
+C triangular part of the array is not referenced.
+C On exit, the leading M-by-M upper triangular part (if
+C UPLO = 'U'), or lower triangular part (if UPLO = 'L') of
+C this array contains the corresponding triangular part of
+C _
+C the computed matrix R.
+C
+C LDR INTEGER
+C The leading dimension of array R. LDR >= MAX(1,M).
+C
+C A (input) DOUBLE PRECISION array, dimension (LDA,k), where
+C k = N when SIDE = 'L', and TRANS = 'N', or
+C SIDE = 'R', and TRANS = 'T';
+C k = M when SIDE = 'R', and TRANS = 'N', or
+C SIDE = 'L', and TRANS = 'T'.
+C On entry, if SIDE = 'L', and TRANS = 'N', or
+C SIDE = 'R', and TRANS = 'T',
+C the leading M-by-N part of this array must contain the
+C matrix A.
+C On entry, if SIDE = 'R', and TRANS = 'N', or
+C SIDE = 'L', and TRANS = 'T',
+C the leading N-by-M part of this array must contain the
+C matrix A.
+C
+C LDA INTEGER
+C The leading dimension of array A. LDA >= MAX(1,l), where
+C l = M when SIDE = 'L', and TRANS = 'N', or
+C SIDE = 'R', and TRANS = 'T';
+C l = N when SIDE = 'R', and TRANS = 'N', or
+C SIDE = 'L', and TRANS = 'T'.
+C
+C B (input) DOUBLE PRECISION array, dimension (LDB,p), where
+C p = M when SIDE = 'L';
+C p = N when SIDE = 'R'.
+C On entry, the leading N-by-M part, if SIDE = 'L', or
+C M-by-N part, if SIDE = 'R', of this array must contain the
+C matrix B.
+C
+C LDB INTEGER
+C The leading dimension of array B.
+C LDB >= MAX(1,N), if SIDE = 'L';
+C LDB >= MAX(1,M), if SIDE = 'R'.
+C
+C Error Indicator
+C
+C INFO INTEGER
+C = 0: successful exit;
+C < 0: if INFO = -i, the i-th argument had an illegal
+C value.
+C
+C METHOD
+C
+C The matrix expression is evaluated taking the triangular
+C structure into account. BLAS 2 operations are used. A block
+C algorithm can be easily constructed; it can use BLAS 3 GEMM
+C operations for most computations, and calls of this BLAS 2
+C algorithm for computing the triangles.
+C
+C FURTHER COMMENTS
+C
+C The main application of this routine is when the result should
+C be a symmetric matrix, e.g., when B = X*op( A )', for (1), or
+C B = op( A )'*X, for (2), where B is already available and X = X'.
+C
+C CONTRIBUTORS
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1999.
+C
+C REVISIONS
+C
+C V. Sima, Research Institute for Informatics, Bucharest, Mar. 2004.
+C
+C KEYWORDS
+C
+C Elementary matrix operations, matrix algebra, matrix operations.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+C .. Scalar Arguments ..
+ CHARACTER SIDE, TRANS, UPLO
+ INTEGER INFO, LDA, LDB, LDR, M, N
+ DOUBLE PRECISION ALPHA, BETA
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*), B(LDB,*), R(LDR,*)
+C .. Local Scalars ..
+ LOGICAL LSIDE, LTRANS, LUPLO
+ INTEGER J
+C .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C .. External Subroutines ..
+ EXTERNAL DGEMV, DLASCL, DLASET, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC MAX
+C .. Executable Statements ..
+C
+C Test the input scalar arguments.
+C
+ INFO = 0
+ LSIDE = LSAME( SIDE, 'L' )
+ LUPLO = LSAME( UPLO, 'U' )
+ LTRANS = LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' )
+C
+ IF( ( .NOT.LSIDE ).AND.( .NOT.LSAME( SIDE, 'R' ) ) )THEN
+ INFO = -1
+ ELSE IF( ( .NOT.LUPLO ).AND.( .NOT.LSAME( UPLO, 'L' ) ) )THEN
+ INFO = -2
+ ELSE IF( ( .NOT.LTRANS ).AND.( .NOT.LSAME( TRANS, 'N' ) ) )THEN
+ INFO = -3
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -5
+ ELSE IF( LDR.LT.MAX( 1, M ) ) THEN
+ INFO = -9
+ ELSE IF( LDA.LT.1 .OR.
+ $ ( ( ( LSIDE .AND. .NOT.LTRANS ) .OR.
+ $ ( .NOT.LSIDE .AND. LTRANS ) ) .AND. LDA.LT.M ) .OR.
+ $ ( ( ( LSIDE .AND. LTRANS ) .OR.
+ $ ( .NOT.LSIDE .AND. .NOT.LTRANS ) ) .AND. LDA.LT.N ) ) THEN
+ INFO = -11
+ ELSE IF( LDB.LT.1 .OR.
+ $ ( LSIDE .AND. LDB.LT.N ) .OR.
+ $ ( .NOT.LSIDE .AND. LDB.LT.M ) ) THEN
+ INFO = -13
+ END IF
+C
+ IF ( INFO.NE.0 ) THEN
+C
+C Error return.
+C
+ CALL XERBLA( 'MB01RX', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ IF ( M.EQ.0 )
+ $ RETURN
+C
+ IF ( BETA.EQ.ZERO .OR. N.EQ.0 ) THEN
+ IF ( ALPHA.EQ.ZERO ) THEN
+C
+C Special case alpha = 0.
+C
+ CALL DLASET( UPLO, M, M, ZERO, ZERO, R, LDR )
+ ELSE
+C
+C Special case beta = 0 or N = 0.
+C
+ IF ( ALPHA.NE.ONE )
+ $ CALL DLASCL( UPLO, 0, 0, ONE, ALPHA, M, M, R, LDR, INFO )
+ END IF
+ RETURN
+ END IF
+C
+C General case: beta <> 0.
+C Compute the required triangle of (1) or (2) using BLAS 2
+C operations.
+C
+ IF( LSIDE ) THEN
+ IF( LUPLO ) THEN
+ IF ( LTRANS ) THEN
+ DO 10 J = 1, M
+ CALL DGEMV( TRANS, N, J, BETA, A, LDA, B(1,J), 1,
+ $ ALPHA, R(1,J), 1 )
+ 10 CONTINUE
+ ELSE
+ DO 20 J = 1, M
+ CALL DGEMV( TRANS, J, N, BETA, A, LDA, B(1,J), 1,
+ $ ALPHA, R(1,J), 1 )
+ 20 CONTINUE
+ END IF
+ ELSE
+ IF ( LTRANS ) THEN
+ DO 30 J = 1, M
+ CALL DGEMV( TRANS, N, M-J+1, BETA, A(1,J), LDA,
+ $ B(1,J), 1, ALPHA, R(J,J), 1 )
+ 30 CONTINUE
+ ELSE
+ DO 40 J = 1, M
+ CALL DGEMV( TRANS, M-J+1, N, BETA, A(J,1), LDA,
+ $ B(1,J), 1, ALPHA, R(J,J), 1 )
+ 40 CONTINUE
+ END IF
+ END IF
+C
+ ELSE
+ IF( LUPLO ) THEN
+ IF( LTRANS ) THEN
+ DO 50 J = 1, M
+ CALL DGEMV( 'NoTranspose', J, N, BETA, B, LDB, A(J,1),
+ $ LDA, ALPHA, R(1,J), 1 )
+ 50 CONTINUE
+ ELSE
+ DO 60 J = 1, M
+ CALL DGEMV( 'NoTranspose', J, N, BETA, B, LDB, A(1,J),
+ $ 1, ALPHA, R(1,J), 1 )
+ 60 CONTINUE
+ END IF
+ ELSE
+ IF( LTRANS ) THEN
+ DO 70 J = 1, M
+ CALL DGEMV( 'NoTranspose', M-J+1, N, BETA, B(J,1),
+ $ LDB, A(J,1), LDA, ALPHA, R(J,J), 1 )
+ 70 CONTINUE
+ ELSE
+ DO 80 J = 1, M
+ CALL DGEMV( 'NoTranspose', M-J+1, N, BETA, B(J,1),
+ $ LDB, A(1,J), 1, ALPHA, R(J,J), 1 )
+ 80 CONTINUE
+ END IF
+ END IF
+ END IF
+C
+ RETURN
+C *** Last line of MB01RX ***
+ END
Property changes on: trunk/octave-forge/main/control/devel/ncfsyn/MB01RX.f
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/control/devel/ncfsyn/MB01RY.f
===================================================================
--- trunk/octave-forge/main/control/devel/ncfsyn/MB01RY.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ncfsyn/MB01RY.f 2011-07-31 19:51:16 UTC (rev 8425)
@@ -0,0 +1,429 @@
+ SUBROUTINE MB01RY( SIDE, UPLO, TRANS, M, ALPHA, BETA, R, LDR, H,
+ $ LDH, B, LDB, DWORK, INFO )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2010 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To compute either the upper or lower triangular part of one of the
+C matrix formulas
+C _
+C R = alpha*R + beta*op( H )*B, (1)
+C _
+C R = alpha*R + beta*B*op( H ), (2)
+C _
+C where alpha and beta are scalars, H, B, R, and R are m-by-m
+C matrices, H is an upper Hessenberg matrix, and op( H ) is one of
+C
+C op( H ) = H or op( H ) = H', the transpose of H.
+C
+C The result is overwritten on R.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C SIDE CHARACTER*1
+C Specifies whether the Hessenberg matrix H appears on the
+C left or right in the matrix product as follows:
+C _
+C = 'L': R = alpha*R + beta*op( H )*B;
+C _
+C = 'R': R = alpha*R + beta*B*op( H ).
+C
+C UPLO CHARACTER*1 _
+C Specifies which triangles of the matrices R and R are
+C computed and given, respectively, as follows:
+C = 'U': the upper triangular part;
+C = 'L': the lower triangular part.
+C
+C TRANS CHARACTER*1
+C Specifies the form of op( H ) to be used in the matrix
+C multiplication as follows:
+C = 'N': op( H ) = H;
+C = 'T': op( H ) = H';
+C = 'C': op( H ) = H'.
+C
+C Input/Output Parameters
+C
+C M (input) INTEGER _
+C The order of the matrices R, R, H and B. M >= 0.
+C
+C ALPHA (input) DOUBLE PRECISION
+C The scalar alpha. When alpha is zero then R need not be
+C set before entry.
+C
+C BETA (input) DOUBLE PRECISION
+C The scalar beta. When beta is zero then H and B are not
+C referenced.
+C
+C R (input/output) DOUBLE PRECISION array, dimension (LDR,M)
+C On entry with UPLO = 'U', the leading M-by-M upper
+C triangular part of this array must contain the upper
+C triangular part of the matrix R; the strictly lower
+C triangular part of the array is not referenced.
+C On entry with UPLO = 'L', the leading M-by-M lower
+C triangular part of this array must contain the lower
+C triangular part of the matrix R; the strictly upper
+C triangular part of the array is not referenced.
+C On exit, the leading M-by-M upper triangular part (if
+C UPLO = 'U'), or lower triangular part (if UPLO = 'L') of
+C this array contains the corresponding triangular part of
+C _
+C the computed matrix R.
+C
+C LDR INTEGER
+C The leading dimension of array R. LDR >= MAX(1,M).
+C
+C H (input) DOUBLE PRECISION array, dimension (LDH,M)
+C On entry, the leading M-by-M upper Hessenberg part of
+C this array must contain the upper Hessenberg part of the
+C matrix H.
+C The elements below the subdiagonal are not referenced,
+C except possibly for those in the first column, which
+C could be overwritten, but are restored on exit.
+C
+C LDH INTEGER
+C The leading dimension of array H. LDH >= MAX(1,M).
+C
+C B (input) DOUBLE PRECISION array, dimension (LDB,M)
+C On entry, the leading M-by-M part of this array must
+C contain the matrix B.
+C
+C LDB INTEGER
+C The leading dimension of array B. LDB >= MAX(1,M).
+C
+C Workspace
+C
+C DWORK DOUBLE PRECISION array, dimension (LDWORK)
+C LDWORK >= M, if beta <> 0 and SIDE = 'L';
+C LDWORK >= 0, if beta = 0 or SIDE = 'R'.
+C This array is not referenced when beta = 0 or SIDE = 'R'.
+C
+C Error Indicator
+C
+C INFO INTEGER
+C = 0: successful exit;
+C < 0: if INFO = -i, the i-th argument had an illegal
+C value.
+C
+C METHOD
+C
+C The matrix expression is efficiently evaluated taking the
+C Hessenberg/triangular structure into account. BLAS 2 operations
+C are used. A block algorithm can be constructed; it can use BLAS 3
+C GEMM operations for most computations, and calls of this BLAS 2
+C algorithm for computing the triangles.
+C
+C FURTHER COMMENTS
+C
+C The main application of this routine is when the result should
+C be a symmetric matrix, e.g., when B = X*op( H )', for (1), or
+C B = op( H )'*X, for (2), where B is already available and X = X'.
+C
+C CONTRIBUTORS
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1999.
+C
+C REVISIONS
+C
+C -
+C
+C KEYWORDS
+C
+C Elementary matrix operations, matrix algebra, matrix operations.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+C .. Scalar Arguments ..
+ CHARACTER SIDE, TRANS, UPLO
+ INTEGER INFO, LDB, LDH, LDR, M
+ DOUBLE PRECISION ALPHA, BETA
+C .. Array Arguments ..
+ DOUBLE PRECISION B(LDB,*), DWORK(*), H(LDH,*), R(LDR,*)
+C .. Local Scalars ..
+ LOGICAL LSIDE, LTRANS, LUPLO
+ INTEGER I, J
+C .. External Functions ..
+ LOGICAL LSAME
+ DOUBLE PRECISION DDOT
+ EXTERNAL DDOT, LSAME
+C .. External Subroutines ..
+ EXTERNAL DCOPY, DGEMV, DLASCL, DLASET, DSCAL, DSWAP,
+ $ DTRMV, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC MAX, MIN
+C .. Executable Statements ..
+C
+C Test the input scalar arguments.
+C
+ INFO = 0
+ LSIDE = LSAME( SIDE, 'L' )
+ LUPLO = LSAME( UPLO, 'U' )
+ LTRANS = LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' )
+C
+ IF( ( .NOT.LSIDE ).AND.( .NOT.LSAME( SIDE, 'R' ) ) )THEN
+ INFO = -1
+ ELSE IF( ( .NOT.LUPLO ).AND.( .NOT.LSAME( UPLO, 'L' ) ) )THEN
+ INFO = -2
+ ELSE IF( ( .NOT.LTRANS ).AND.( .NOT.LSAME( TRANS, 'N' ) ) )THEN
+ INFO = -3
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDR.LT.MAX( 1, M ) ) THEN
+ INFO = -8
+ ELSE IF( LDH.LT.MAX( 1, M ) ) THEN
+ INFO = -10
+ ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
+ INFO = -12
+ END IF
+C
+ IF ( INFO.NE.0 ) THEN
+C
+C Error return.
+C
+ CALL XERBLA( 'MB01RY', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ IF ( M.EQ.0 )
+ $ RETURN
+C
+ IF ( BETA.EQ.ZERO ) THEN
+ IF ( ALPHA.EQ.ZERO ) THEN
+C
+C Special case when both alpha = 0 and beta = 0.
+C
+ CALL DLASET( UPLO, M, M, ZERO, ZERO, R, LDR )
+ ELSE
+C
+C Special case beta = 0.
+C
+ IF ( ALPHA.NE.ONE )
+ $ CALL DLASCL( UPLO, 0, 0, ONE, ALPHA, M, M, R, LDR, INFO )
+ END IF
+ RETURN
+ END IF
+C
+C General case: beta <> 0.
+C Compute the required triangle of (1) or (2) using BLAS 2
+C operations.
+C
+ IF( LSIDE ) THEN
+C
+C To avoid repeated references to the subdiagonal elements of H,
+C these are swapped with the corresponding elements of H in the
+C first column, and are finally restored.
+C
+ IF( M.GT.2 )
+ $ CALL DSWAP( M-2, H( 3, 2 ), LDH+1, H( 3, 1 ), 1 )
+C
+ IF( LUPLO ) THEN
+ IF ( LTRANS ) THEN
+C
+ DO 20 J = 1, M
+C
+C Multiply the transposed upper triangle of the leading
+C j-by-j submatrix of H by the leading part of the j-th
+C column of B.
+C
+ CALL DCOPY( J, B( 1, J ), 1, DWORK, 1 )
+ CALL DTRMV( 'Upper', TRANS, 'Non-unit', J, H, LDH,
+ $ DWORK, 1 )
+C
+C Add the contribution of the subdiagonal of H to
+C the j-th column of the product.
+C
+ DO 10 I = 1, MIN( J, M - 1 )
+ R( I, J ) = ALPHA*R( I, J ) + BETA*( DWORK( I ) +
+ $ H( I+1, 1 )*B( I+1, J ) )
+ 10 CONTINUE
+C
+ 20 CONTINUE
+C
+ R( M, M ) = ALPHA*R( M, M ) + BETA*DWORK( M )
+C
+ ELSE
+C
+ DO 40 J = 1, M
+C
+C Multiply the upper triangle of the leading j-by-j
+C submatrix of H by the leading part of the j-th column
+C of B.
+C
+ CALL DCOPY( J, B( 1, J ), 1, DWORK, 1 )
+ CALL DTRMV( 'Upper', TRANS, 'Non-unit', J, H, LDH,
+ $ DWORK, 1 )
+ IF( J.LT.M ) THEN
+C
+C Multiply the remaining right part of the leading
+C j-by-M submatrix of H by the trailing part of the
+C j-th column of B.
+C
+ CALL DGEMV( TRANS, J, M-J, BETA, H( 1, J+1 ), LDH,
+ $ B( J+1, J ), 1, ALPHA, R( 1, J ), 1 )
+ ELSE
+ CALL DSCAL( M, ALPHA, R( 1, M ), 1 )
+ END IF
+C
+C Add the contribution of the subdiagonal of H to
+C the j-th column of the product.
+C
+ R( 1, J ) = R( 1, J ) + BETA*DWORK( 1 )
+C
+ DO 30 I = 2, J
+ R( I, J ) = R( I, J ) + BETA*( DWORK( I ) +
+ $ H( I, 1 )*B( I-1, J ) )
+ 30 CONTINUE
+C
+ 40 CONTINUE
+C
+ END IF
+C
+ ELSE
+C
+ IF ( LTRANS ) THEN
+C
+ DO 60 J = M, 1, -1
+C
+C Multiply the transposed upper triangle of the trailing
+C (M-j+1)-by-(M-j+1) submatrix of H by the trailing part
+C of the j-th column of B.
+C
+ CALL DCOPY( M-J+1, B( J, J ), 1, DWORK( J ), 1 )
+ CALL DTRMV( 'Upper', TRANS, 'Non-unit', M-J+1,
+ $ H( J, J ), LDH, DWORK( J ), 1 )
+ IF( J.GT.1 ) THEN
+C
+C Multiply the remaining left part of the trailing
+C (M-j+1)-by-(j-1) submatrix of H' by the leading
+C part of the j-th column of B.
+C
+ CALL DGEMV( TRANS, J-1, M-J+1, BETA, H( 1, J ),
+ $ LDH, B( 1, J ), 1, ALPHA, R( J, J ),
+ $ 1 )
+ ELSE
+ CALL DSCAL( M, ALPHA, R( 1, 1 ), 1 )
+ END IF
+C
+C Add the contribution of the subdiagonal of H to
+C the j-th column of the product.
+C
+ DO 50 I = J, M - 1
+ R( I, J ) = R( I, J ) + BETA*( DWORK( I ) +
+ $ H( I+1, 1 )*B( I+1, J ) )
+ 50 CONTINUE
+C
+ R( M, J ) = R( M, J ) + BETA*DWORK( M )
+ 60 CONTINUE
+C
+ ELSE
+C
+ DO 80 J = M, 1, -1
+C
+C Multiply the upper triangle of the trailing
+C (M-j+1)-by-(M-j+1) submatrix of H by the trailing
+C part of the j-th column of B.
+C
+ CALL DCOPY( M-J+1, B( J, J ), 1, DWORK( J ), 1 )
+ CALL DTRMV( 'Upper', TRANS, 'Non-unit', M-J+1,
+ $ H( J, J ), LDH, DWORK( J ), 1 )
+C
+C Add the contribution of the subdiagonal of H to
+C the j-th column of the product.
+C
+ DO 70 I = MAX( J, 2 ), M
+ R( I, J ) = ALPHA*R( I, J ) + BETA*( DWORK( I )
+ $ + H( I, 1 )*B( I-1, J ) )
+ 70 CONTINUE
+C
+ 80 CONTINUE
+C
+ R( 1, 1 ) = ALPHA*R( 1, 1 ) + BETA*DWORK( 1 )
+C
+ END IF
+ END IF
+C
+ IF( M.GT.2 )
+ $ CALL DSWAP( M-2, H( 3, 2 ), LDH+1, H( 3, 1 ), 1 )
+C
+ ELSE
+C
+C Row-wise calculations are used for H, if SIDE = 'R' and
+C TRANS = 'T'.
+C
+ IF( LUPLO ) THEN
+ IF( LTRANS ) THEN
+ R( 1, 1 ) = ALPHA*R( 1, 1 ) +
+ $ BETA*DDOT( M, B, LDB, H, LDH )
+C
+ DO 90 J = 2, M
+ CALL DGEMV( 'NoTranspose', J, M-J+2, BETA,
+ $ B( 1, J-1 ), LDB, H( J, J-1 ), LDH,
+ $ ALPHA, R( 1, J ), 1 )
+ 90 CONTINUE
+C
+ ELSE
+C
+ DO 100 J = 1, M - 1
+ CALL DGEMV( 'NoTranspose', J, J+1, BETA, B, LDB,
+ $ H( 1, J ), 1, ALPHA, R( 1, J ), 1 )
+ 100 CONTINUE
+C
+ CALL DGEMV( 'NoTranspose', M, M, BETA, B, LDB,
+ $ H( 1, M ), 1, ALPHA, R( 1, M ), 1 )
+C
+ END IF
+C
+ ELSE
+C
+ IF( LTRANS ) THEN
+C
+ CALL DGEMV( 'NoTranspose', M, M, BETA, B, LDB, H, LDH,
+ $ ALPHA, R( 1, 1 ), 1 )
+C
+ DO 110 J = 2, M
+ CALL DGEMV( 'NoTranspose', M-J+1, M-J+2, BETA,
+ $ B( J, J-1 ), LDB, H( J, J-1 ), LDH, ALPHA,
+ $ R( J, J ), 1 )
+ 110 CONTINUE
+C
+ ELSE
+C
+ DO 120 J = 1, M - 1
+ CALL DGEMV( 'NoTranspose', M-J+1, J+1, BETA,
+ $ B( J, 1 ), LDB, H( 1, J ), 1, ALPHA,
+ $ R( J, J ), 1 )
+ 120 CONTINUE
+C
+ R( M, M ) = ALPHA*R( M, M ) +
+ $ BETA*DDOT( M, B( M, 1 ), LDB, H( 1, M ), 1 )
+C
+ END IF
+ END IF
+ END IF
+C
+ RETURN
+C *** Last line of MB01RY ***
+ END
Property changes on: trunk/octave-forge/main/control/devel/ncfsyn/MB01RY.f
___________________________________________________________________
Added: svn:executable
+ *
Added: trunk/octave-forge/main/control/devel/ncfsyn/MB01SD.f
===================================================================
--- trunk/octave-forge/main/control/devel/ncfsyn/MB01SD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ncfsyn/MB01SD.f 2011-07-31 19:51:16 UTC (rev 8425)
@@ -0,0 +1,123 @@
+ SUBROUTINE MB01SD( JOBS, M, N, A, LDA, R, C )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2010 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To scale a general M-by-N matrix A using the row and column
+C scaling factors in the vectors R and C.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C JOBS CHARACTER*1
+C Specifies the scaling operation to be done, as follows:
+C = 'R': row scaling, i.e., A will be premultiplied
+C by diag(R);
+C = 'C': column scaling, i.e., A will be postmultiplied
+C by diag(C);
+C = 'B': both row and column scaling, i.e., A will be
+C replaced by diag(R) * A * diag(C).
+C
+C Input/Output Parameters
+C
+C M (input) INTEGER
+C The number of rows of the matrix A. M >= 0.
+C
+C N (input) INTEGER
+C The number of columns of the matrix A. N >= 0.
+C
+C A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+C On entry, the M-by-N matrix A.
+C On exit, the scaled matrix. See JOBS for the form of the
+C scaled matrix.
+C
+C LDA INTEGER
+C The leading dimension of the array A. LDA >= max(1,M).
+C
+C R (input) DOUBLE PRECISION array, dimension (M)
+C The row scale factors for A.
+C R is not referenced if JOBS = 'C'.
+C
+C C (input) DOUBLE PRECISION array, dimension (N)
+C The column scale factors for A.
+C C is not referenced if JOBS = 'R'.
+C
+C
+C CONTRIBUTOR
+C
+C A. Varga, German Aerospace Center,
+C DLR Oberpfaffenhofen, April 1998.
+C Based on the RASP routine DMSCAL.
+C
+C ******************************************************************
+C
+C .. Scalar Arguments ..
+ CHARACTER JOBS
+ INTEGER LDA, M, N
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*), C(*), R(*)
+C .. Local Scalars ..
+ INTEGER I, J
+ DOUBLE PRECISION CJ
+C .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C .. Executable Statements ..
+C
+C Quick return if possible.
+C
+ IF( M.EQ.0 .OR. N.EQ.0 )
+ $ RETURN
+C
+ IF( LSAME( JOBS, 'C' ) ) ...
[truncated message content] |
