crsf2csf.f

Go to the documentation of this file.
00001 c Copyright (C) 2010-2012  VZLU Prague, a.s., Czech Republic
00002 c
00003 c Author: Jaroslav Hajek <highegg@gmail.com>
00004 c
00005 c This file is part of Octave.
00006 c
00007 c Octave is free software; you can redistribute it and/or modify
00008 c it under the terms of the GNU General Public License as published by
00009 c the Free Software Foundation; either version 3 of the License, or
00010 c (at your option) any later version.
00011 c
00012 c This program is distributed in the hope that it will be useful,
00013 c but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 c GNU General Public License for more details.
00016 c
00017 c You should have received a copy of the GNU General Public License
00018 c along with this software; see the file COPYING.  If not, see
00019 c <http://www.gnu.org/licenses/>.
00020 c
00021 
00022        subroutine crsf2csf(n,t,u,c,s)
00023        integer n
00024        complex t(n,n),u(n,n)
00025        real c(n-1),s(n-1)
00026        real x,y,z
00027        integer j
00028        do j = 1,n-1
00029           c(j) = 1
00030        end do
00031        j = 1
00032        do while (j < n)
00033 c apply previous rotations to rows
00034          call crcrot1(j,t(1,j),c,s)
00035 
00036          y = t(j+1,j)
00037          if (y /= 0) then
00038 c 2x2 block, form Givens rotation [c, i*s; i*s, c]
00039            z = t(j,j+1)
00040            c(j) = sqrt(z/(z-y))
00041            s(j) = sqrt(y/(y-z))
00042 c apply new rotation to t(j:j+1,j)
00043            call crcrot1(2,t(j,j),c(j),s(j))
00044 c apply all rotations to t(1:j+1,j+1)
00045            call crcrot1(j+1,t(1,j+1),c,s)
00046 c apply new rotation to columns j,j+1
00047            call crcrot2(j+1,t(1,j),t(1,j+1),c(j),s(j))
00048 c zero subdiagonal entry, skip next row
00049            t(j+1,j) = 0
00050            j = j + 2
00051          else
00052            j = j + 1
00053          end if
00054        end do
00055 
00056 c apply rotations to last column if needed
00057        if (j == n) then
00058          call crcrot1(j,t(1,j),c,s)
00059        end if
00060 
00061 c apply stored rotations to all columns of u
00062        do j = 1,n-1
00063          if (c(j) /= 1) then
00064            call crcrot2(n,u(1,j),u(1,j+1),c(j),s(j))
00065          end if
00066        end do
00067 
00068        end subroutine
00069 
00070        subroutine crcrot1(n,x,c,s)
00071 c apply rotations to a column from the left
00072        integer n
00073        complex x(n), t
00074        real c(n-1),s(n-1)
00075        integer i
00076        do i = 1,n-1
00077          if (c(i) /= 1) then
00078            t = x(i)*c(i) - x(i+1)*cmplx(0,s(i))
00079            x(i+1) = x(i+1)*c(i) - x(i)*cmplx(0,s(i))
00080            x(i) = t
00081          endif
00082        end do
00083        end subroutine
00084 
00085        subroutine crcrot2(n,x,y,c,s)
00086 c apply a single rotation from the right to a pair of columns
00087        integer n
00088        complex x(n),y(n),t
00089        real c, s
00090        integer i
00091        do i = 1,n
00092          t = x(i)*c + y(i)*cmplx(0,s)
00093          y(i) = y(i)*c + x(i)*cmplx(0,s)
00094          x(i) = t
00095        end do
00096        end subroutine
00097 
00098 
00099 
00100 
00101 
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines