/************************************************************************* Copyright (c) 1992-2007 The University of Tennessee. All rights reserved. Contributors: * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to pseudocode. See subroutines comments for additional copyrights. >>> SOURCE LICENSE >>> This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation (www.fsf.org); either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A copy of the GNU General Public License is available at http://www.fsf.org/licensing/licenses >>> END OF LICENSE >>> *************************************************************************/ using System; namespace alglib { public class rotations { /************************************************************************* Application of a sequence of elementary rotations to a matrix The algorithm pre-multiplies the matrix by a sequence of rotation transformations which is given by arrays C and S. Depending on the value of the IsForward parameter either 1 and 2, 3 and 4 and so on (if IsForward=true) rows are rotated, or the rows N and N-1, N-2 and N-3 and so on, are rotated. Not the whole matrix but only a part of it is transformed (rows from M1 to M2, columns from N1 to N2). Only the elements of this submatrix are changed. Input parameters: IsForward - the sequence of the rotation application. M1,M2 - the range of rows to be transformed. N1, N2 - the range of columns to be transformed. C,S - transformation coefficients. Array whose index ranges within [1..M2-M1]. A - processed matrix. WORK - working array whose index ranges within [N1..N2]. Output parameters: A - transformed matrix. Utility subroutine. *************************************************************************/ public static void applyrotationsfromtheleft(bool isforward, int m1, int m2, int n1, int n2, ref double[] c, ref double[] s, ref double[,] a, ref double[] work) { int j = 0; int jp1 = 0; double ctemp = 0; double stemp = 0; double temp = 0; int i_ = 0; if( m1>m2 | n1>n2 ) { return; } // // Form P * A // if( isforward ) { if( n1!=n2 ) { // // Common case: N1<>N2 // for(j=m1; j<=m2-1; j++) { ctemp = c[j-m1+1]; stemp = s[j-m1+1]; if( ctemp!=1 | stemp!=0 ) { jp1 = j+1; for(i_=n1; i_<=n2;i_++) { work[i_] = ctemp*a[jp1,i_]; } for(i_=n1; i_<=n2;i_++) { work[i_] = work[i_] - stemp*a[j,i_]; } for(i_=n1; i_<=n2;i_++) { a[j,i_] = ctemp*a[j,i_]; } for(i_=n1; i_<=n2;i_++) { a[j,i_] = a[j,i_] + stemp*a[jp1,i_]; } for(i_=n1; i_<=n2;i_++) { a[jp1,i_] = work[i_]; } } } } else { // // Special case: N1=N2 // for(j=m1; j<=m2-1; j++) { ctemp = c[j-m1+1]; stemp = s[j-m1+1]; if( ctemp!=1 | stemp!=0 ) { temp = a[j+1,n1]; a[j+1,n1] = ctemp*temp-stemp*a[j,n1]; a[j,n1] = stemp*temp+ctemp*a[j,n1]; } } } } else { if( n1!=n2 ) { // // Common case: N1<>N2 // for(j=m2-1; j>=m1; j--) { ctemp = c[j-m1+1]; stemp = s[j-m1+1]; if( ctemp!=1 | stemp!=0 ) { jp1 = j+1; for(i_=n1; i_<=n2;i_++) { work[i_] = ctemp*a[jp1,i_]; } for(i_=n1; i_<=n2;i_++) { work[i_] = work[i_] - stemp*a[j,i_]; } for(i_=n1; i_<=n2;i_++) { a[j,i_] = ctemp*a[j,i_]; } for(i_=n1; i_<=n2;i_++) { a[j,i_] = a[j,i_] + stemp*a[jp1,i_]; } for(i_=n1; i_<=n2;i_++) { a[jp1,i_] = work[i_]; } } } } else { // // Special case: N1=N2 // for(j=m2-1; j>=m1; j--) { ctemp = c[j-m1+1]; stemp = s[j-m1+1]; if( ctemp!=1 | stemp!=0 ) { temp = a[j+1,n1]; a[j+1,n1] = ctemp*temp-stemp*a[j,n1]; a[j,n1] = stemp*temp+ctemp*a[j,n1]; } } } } } /************************************************************************* Application of a sequence of elementary rotations to a matrix The algorithm post-multiplies the matrix by a sequence of rotation transformations which is given by arrays C and S. Depending on the value of the IsForward parameter either 1 and 2, 3 and 4 and so on (if IsForward=true) rows are rotated, or the rows N and N-1, N-2 and N-3 and so on are rotated. Not the whole matrix but only a part of it is transformed (rows from M1 to M2, columns from N1 to N2). Only the elements of this submatrix are changed. Input parameters: IsForward - the sequence of the rotation application. M1,M2 - the range of rows to be transformed. N1, N2 - the range of columns to be transformed. C,S - transformation coefficients. Array whose index ranges within [1..N2-N1]. A - processed matrix. WORK - working array whose index ranges within [M1..M2]. Output parameters: A - transformed matrix. Utility subroutine. *************************************************************************/ public static void applyrotationsfromtheright(bool isforward, int m1, int m2, int n1, int n2, ref double[] c, ref double[] s, ref double[,] a, ref double[] work) { int j = 0; int jp1 = 0; double ctemp = 0; double stemp = 0; double temp = 0; int i_ = 0; // // Form A * P' // if( isforward ) { if( m1!=m2 ) { // // Common case: M1<>M2 // for(j=n1; j<=n2-1; j++) { ctemp = c[j-n1+1]; stemp = s[j-n1+1]; if( ctemp!=1 | stemp!=0 ) { jp1 = j+1; for(i_=m1; i_<=m2;i_++) { work[i_] = ctemp*a[i_,jp1]; } for(i_=m1; i_<=m2;i_++) { work[i_] = work[i_] - stemp*a[i_,j]; } for(i_=m1; i_<=m2;i_++) { a[i_,j] = ctemp*a[i_,j]; } for(i_=m1; i_<=m2;i_++) { a[i_,j] = a[i_,j] + stemp*a[i_,jp1]; } for(i_=m1; i_<=m2;i_++) { a[i_,jp1] = work[i_]; } } } } else { // // Special case: M1=M2 // for(j=n1; j<=n2-1; j++) { ctemp = c[j-n1+1]; stemp = s[j-n1+1]; if( ctemp!=1 | stemp!=0 ) { temp = a[m1,j+1]; a[m1,j+1] = ctemp*temp-stemp*a[m1,j]; a[m1,j] = stemp*temp+ctemp*a[m1,j]; } } } } else { if( m1!=m2 ) { // // Common case: M1<>M2 // for(j=n2-1; j>=n1; j--) { ctemp = c[j-n1+1]; stemp = s[j-n1+1]; if( ctemp!=1 | stemp!=0 ) { jp1 = j+1; for(i_=m1; i_<=m2;i_++) { work[i_] = ctemp*a[i_,jp1]; } for(i_=m1; i_<=m2;i_++) { work[i_] = work[i_] - stemp*a[i_,j]; } for(i_=m1; i_<=m2;i_++) { a[i_,j] = ctemp*a[i_,j]; } for(i_=m1; i_<=m2;i_++) { a[i_,j] = a[i_,j] + stemp*a[i_,jp1]; } for(i_=m1; i_<=m2;i_++) { a[i_,jp1] = work[i_]; } } } } else { // // Special case: M1=M2 // for(j=n2-1; j>=n1; j--) { ctemp = c[j-n1+1]; stemp = s[j-n1+1]; if( ctemp!=1 | stemp!=0 ) { temp = a[m1,j+1]; a[m1,j+1] = ctemp*temp-stemp*a[m1,j]; a[m1,j] = stemp*temp+ctemp*a[m1,j]; } } } } } /************************************************************************* The subroutine generates the elementary rotation, so that: [ CS SN ] . [ F ] = [ R ] [ -SN CS ] [ G ] [ 0 ] CS**2 + SN**2 = 1 *************************************************************************/ public static void generaterotation(double f, double g, ref double cs, ref double sn, ref double r) { double f1 = 0; double g1 = 0; if( g==0 ) { cs = 1; sn = 0; r = f; } else { if( f==0 ) { cs = 0; sn = 1; r = g; } else { f1 = f; g1 = g; r = Math.Sqrt(AP.Math.Sqr(f1)+AP.Math.Sqr(g1)); cs = f1/r; sn = g1/r; if( Math.Abs(f)>Math.Abs(g) & cs<0 ) { cs = -cs; sn = -sn; r = -r; } } } } private static void testrotations() { double[,] al1 = new double[0,0]; double[,] al2 = new double[0,0]; double[,] ar1 = new double[0,0]; double[,] ar2 = new double[0,0]; double[] cl = new double[0]; double[] sl = new double[0]; double[] cr = new double[0]; double[] sr = new double[0]; double[] w = new double[0]; int m = 0; int n = 0; int maxmn = 0; double t = 0; int pass = 0; int passcount = 0; int i = 0; int j = 0; double err = 0; double maxerr = 0; bool isforward = new bool(); passcount = 1000; maxerr = 0; for(pass=1; pass<=passcount; pass++) { // // settings // m = 2+AP.Math.RandomInteger(50); n = 2+AP.Math.RandomInteger(50); isforward = AP.Math.RandomReal()>0.5; maxmn = Math.Max(m, n); al1 = new double[m+1, n+1]; al2 = new double[m+1, n+1]; ar1 = new double[m+1, n+1]; ar2 = new double[m+1, n+1]; cl = new double[m-1+1]; sl = new double[m-1+1]; cr = new double[n-1+1]; sr = new double[n-1+1]; w = new double[maxmn+1]; // // matrices and rotaions // for(i=1; i<=m; i++) { for(j=1; j<=n; j++) { al1[i,j] = 2*AP.Math.RandomReal()-1; al2[i,j] = al1[i,j]; ar1[i,j] = al1[i,j]; ar2[i,j] = al1[i,j]; } } for(i=1; i<=m-1; i++) { t = 2*Math.PI*AP.Math.RandomReal(); cl[i] = Math.Cos(t); sl[i] = Math.Sin(t); } for(j=1; j<=n-1; j++) { t = 2*Math.PI*AP.Math.RandomReal(); cr[j] = Math.Cos(t); sr[j] = Math.Sin(t); } // // Test left // applyrotationsfromtheleft(isforward, 1, m, 1, n, ref cl, ref sl, ref al1, ref w); for(j=1; j<=n; j++) { applyrotationsfromtheleft(isforward, 1, m, j, j, ref cl, ref sl, ref al2, ref w); } err = 0; for(i=1; i<=m; i++) { for(j=1; j<=n; j++) { err = Math.Max(err, Math.Abs(al1[i,j]-al2[i,j])); } } maxerr = Math.Max(err, maxerr); // // Test right // applyrotationsfromtheright(isforward, 1, m, 1, n, ref cr, ref sr, ref ar1, ref w); for(i=1; i<=m; i++) { applyrotationsfromtheright(isforward, i, i, 1, n, ref cr, ref sr, ref ar2, ref w); } err = 0; for(i=1; i<=m; i++) { for(j=1; j<=n; j++) { err = Math.Max(err, Math.Abs(ar1[i,j]-ar2[i,j])); } } maxerr = Math.Max(err, maxerr); } System.Console.Write("TESTING ROTATIONS"); System.Console.WriteLine(); System.Console.Write("Pass count "); System.Console.Write("{0,0:d}",passcount); System.Console.WriteLine(); System.Console.Write("Error is "); System.Console.Write("{0,5:E3}",maxerr); System.Console.WriteLine(); } } }