Mercurial > hg > octave-lyh

changeset 7886:e3e94982dfd4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Convert qrshift, qrinsert, qrdelete, qrupdate to allow single precision arguments
author David Bateman <dbateman@free.fr>
date Sun, 08 Jun 2008 20:41:01 +0200
parents f336dd8e96d0
children 627b10572d82
files src/ChangeLog src/DLD-FUNCTIONS/qr.cc
diffstat 2 files changed, 444 insertions(+), 195 deletions(-) [+]
line wrap: on
line diff
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,8 @@
+2008-06-12  David Bateman  <dbateman@free.fr>
+
+	* DLD-FUNCTIONS/qr.cc (Fqrupdate, Fqrinsert, Fqrshift, Fqrdelete):
+	Allow single precision arguments, add tests for single precision.
+
 2008-06-11  John W. Eaton  <jwe@octave.org>
 
 	* ov-base.cc (octave_base_value::streamoff_value,
--- a/src/DLD-FUNCTIONS/qr.cc
+++ b/src/DLD-FUNCTIONS/qr.cc
@@ -773,31 +773,69 @@
 	      && argu.is_real_matrix () 
 	      && argv.is_real_matrix ())
             {
-              // all real case
-              Matrix Q = argq.matrix_value ();
-              Matrix R = argr.matrix_value ();
-              Matrix u = argu.matrix_value ();
-              Matrix v = argv.matrix_value ();
+	      // all real case
+	      if (argq.is_single_type () 
+		  || argr.is_single_type () 
+		  || argu.is_single_type () 
+		  || argv.is_single_type ())
+		{
+		  FloatMatrix Q = argq.float_matrix_value ();
+		  FloatMatrix R = argr.float_matrix_value ();
+		  FloatMatrix u = argu.float_matrix_value ();
+		  FloatMatrix v = argv.float_matrix_value ();
+
+		  FloatQR fact (Q, R);
+		  fact.update (u, v);
 
-              QR fact (Q, R);
-              fact.update (u, v);
+		  retval(1) = fact.R ();
+		  retval(0) = fact.Q ();
+		}
+	      else
+		{
+		  Matrix Q = argq.matrix_value ();
+		  Matrix R = argr.matrix_value ();
+		  Matrix u = argu.matrix_value ();
+		  Matrix v = argv.matrix_value ();
 
-              retval(1) = fact.R ();
-              retval(0) = fact.Q ();
+		  QR fact (Q, R);
+		  fact.update (u, v);
+
+		  retval(1) = fact.R ();
+		  retval(0) = fact.Q ();
+		}
             }
           else
             {
               // complex case
-              ComplexMatrix Q = argq.complex_matrix_value ();
-              ComplexMatrix R = argr.complex_matrix_value ();
-              ComplexMatrix u = argu.complex_matrix_value ();
-              ComplexMatrix v = argv.complex_matrix_value ();
+	      if (argq.is_single_type () 
+		  || argr.is_single_type () 
+		  || argu.is_single_type () 
+		  || argv.is_single_type ())
+		{
+		  FloatComplexMatrix Q = argq.float_complex_matrix_value ();
+		  FloatComplexMatrix R = argr.float_complex_matrix_value ();
+		  FloatComplexMatrix u = argu.float_complex_matrix_value ();
+		  FloatComplexMatrix v = argv.float_complex_matrix_value ();
 
-              ComplexQR fact (Q, R);
-              fact.update (u, v);
+		  FloatComplexQR fact (Q, R);
+		  fact.update (u, v);
               
-              retval(1) = fact.R ();
-              retval(0) = fact.Q ();
+		  retval(1) = fact.R ();
+		  retval(0) = fact.Q ();
+		}
+	      else
+		{
+		  ComplexMatrix Q = argq.complex_matrix_value ();
+		  ComplexMatrix R = argr.complex_matrix_value ();
+		  ComplexMatrix u = argu.complex_matrix_value ();
+		  ComplexMatrix v = argv.complex_matrix_value ();
+
+		  ComplexQR fact (Q, R);
+		  fact.update (u, v);
+              
+		  retval(1) = fact.R ();
+		  retval(0) = fact.Q ();
+		}
             }
         }
       else
@@ -809,7 +847,7 @@
   return retval;
 }
 /*
-%!test
+%!shared A, u, v, Ac, uc, vc
 %! A = [0.091364  0.613038  0.999083;
 %!      0.594638  0.425302  0.603537;
 %!      0.383594  0.291238  0.085574;
@@ -826,6 +864,24 @@
 %!      0.24295;
 %!      0.43167 ];
 %!
+%! Ac = [0.620405 + 0.956953i  0.480013 + 0.048806i  0.402627 + 0.338171i;
+%!      0.589077 + 0.658457i  0.013205 + 0.279323i  0.229284 + 0.721929i;
+%!      0.092758 + 0.345687i  0.928679 + 0.241052i  0.764536 + 0.832406i;
+%!      0.912098 + 0.721024i  0.049018 + 0.269452i  0.730029 + 0.796517i;
+%!      0.112849 + 0.603871i  0.486352 + 0.142337i  0.355646 + 0.151496i ];
+%!
+%! uc = [0.20351 + 0.05401i;
+%!      0.13141 + 0.43708i;
+%!      0.29808 + 0.08789i;
+%!      0.69821 + 0.38844i;
+%!      0.74871 + 0.25821i ];
+%!
+%! vc = [0.85839 + 0.29468i;
+%!      0.20820 + 0.93090i;
+%!      0.86184 + 0.34689i ];
+%!
+
+%!test
 %! [Q,R] = qr(A);
 %! [Q,R] = qrupdate(Q,R,u,v);
 %! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
@@ -833,27 +889,25 @@
 %! assert(norm(vec(Q*R - A - u*v'),Inf) < norm(A)*1e1*eps)
 %! 
 %!test
-%! A = [0.620405 + 0.956953i  0.480013 + 0.048806i  0.402627 + 0.338171i;
-%!      0.589077 + 0.658457i  0.013205 + 0.279323i  0.229284 + 0.721929i;
-%!      0.092758 + 0.345687i  0.928679 + 0.241052i  0.764536 + 0.832406i;
-%!      0.912098 + 0.721024i  0.049018 + 0.269452i  0.730029 + 0.796517i;
-%!      0.112849 + 0.603871i  0.486352 + 0.142337i  0.355646 + 0.151496i ];
-%!
-%! u = [0.20351 + 0.05401i;
-%!      0.13141 + 0.43708i;
-%!      0.29808 + 0.08789i;
-%!      0.69821 + 0.38844i;
-%!      0.74871 + 0.25821i ];
-%!
-%! v = [0.85839 + 0.29468i;
-%!      0.20820 + 0.93090i;
-%!      0.86184 + 0.34689i ];
-%!
-%! [Q,R] = qr(A);
-%! [Q,R] = qrupdate(Q,R,u,v);
+%! [Q,R] = qr(Ac);
+%! [Q,R] = qrupdate(Q,R,uc,vc);
 %! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - A - u*v'),Inf) < norm(A)*1e1*eps)
+%! assert(norm(vec(Q*R - Ac - uc*vc'),Inf) < norm(Ac)*1e1*eps)
+
+%!test
+%! [Q,R] = qr(single(A));
+%! [Q,R] = qrupdate(Q,R,single(u),single(v));
+%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - single(A) - single(u)*single(v)'),Inf) < norm(single(A))*1e1*eps('single'))
+%! 
+%!test
+%! [Q,R] = qr(single(Ac));
+%! [Q,R] = qrupdate(Q,R,single(uc),single(vc));
+%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - single(Ac) - single(uc)*single(vc)'),Inf) < norm(single(Ac))*1e1*eps('single'))
 */
 
 DEFUN_DLD (qrinsert, args, ,
@@ -910,7 +964,7 @@
             && argx.rows () == (row ? 1 : m)
             && argx.columns () == (row ? n : 1))
           {
-            octave_idx_type j = argj.int_value ();
+            octave_idx_type j = argj.idx_type_value ();
 
             if (j >= 1 && j <= (row ? n : m)+1)
               {
@@ -919,36 +973,80 @@
 		    && argx.is_real_matrix ())
                   {
                     // real case
-                    Matrix Q = argq.matrix_value ();
-                    Matrix R = argr.matrix_value ();
-                    Matrix x = argx.matrix_value ();
+		    if (argq.is_single_type () 
+			|| argr.is_single_type () 
+			|| argx.is_single_type ())
+		      {
+			FloatMatrix Q = argq.float_matrix_value ();
+			FloatMatrix R = argr.float_matrix_value ();
+			FloatMatrix x = argx.float_matrix_value ();
 
-                    QR fact (Q, R);
+			FloatQR fact (Q, R);
+
+			if (row) 
+			  fact.insert_row (x, j-1);
+			else 
+			  fact.insert_col (x, j-1);
+
+			retval(1) = fact.R ();
+			retval(0) = fact.Q ();
 
-                    if (row) 
-                      fact.insert_row (x, j-1);
-                    else 
-                      fact.insert_col (x, j-1);
+		      }
+		    else
+		      {
+			Matrix Q = argq.matrix_value ();
+			Matrix R = argr.matrix_value ();
+			Matrix x = argx.matrix_value ();
+
+			QR fact (Q, R);
 
-                    retval(1) = fact.R ();
-                    retval(0) = fact.Q ();
+			if (row) 
+			  fact.insert_row (x, j-1);
+			else 
+			  fact.insert_col (x, j-1);
+
+			retval(1) = fact.R ();
+			retval(0) = fact.Q ();
+
+		      }
                   }
                 else
                   {
                     // complex case
-                    ComplexMatrix Q = argq.complex_matrix_value ();
-                    ComplexMatrix R = argr.complex_matrix_value ();
-                    ComplexMatrix x = argx.complex_matrix_value ();
+		    if (argq.is_single_type () 
+			|| argr.is_single_type () 
+			|| argx.is_single_type ())
+		      {
+			FloatComplexMatrix Q = argq.float_complex_matrix_value ();
+			FloatComplexMatrix R = argr.float_complex_matrix_value ();
+			FloatComplexMatrix x = argx.float_complex_matrix_value ();
 
-                    ComplexQR fact (Q, R);
+			FloatComplexQR fact (Q, R);
+
+			if (row) 
+			  fact.insert_row (x, j-1);
+			else 
+			  fact.insert_col (x, j-1);
 
-                    if (row) 
-                      fact.insert_row (x, j-1);
-                    else 
-                      fact.insert_col (x, j-1);
+			retval(1) = fact.R ();
+			retval(0) = fact.Q ();
+		      }
+		    else
+		      {
+			ComplexMatrix Q = argq.complex_matrix_value ();
+			ComplexMatrix R = argr.complex_matrix_value ();
+			ComplexMatrix x = argx.complex_matrix_value ();
 
-                    retval(1) = fact.R ();
-                    retval(0) = fact.Q ();
+			ComplexQR fact (Q, R);
+
+			if (row) 
+			  fact.insert_row (x, j-1);
+			else 
+			  fact.insert_col (x, j-1);
+
+			retval(1) = fact.R ();
+			retval(0) = fact.Q ();
+		      }
                   }
 
               }
@@ -969,50 +1067,18 @@
 
 /*
 %!test
-%! A = [0.091364  0.613038  0.999083;
-%!      0.594638  0.425302  0.603537;
-%!      0.383594  0.291238  0.085574;
-%!      0.265712  0.268003  0.238409;
-%!      0.669966  0.743851  0.445057 ];
-%!
-%! x = [0.85082;  
-%!      0.76426;  
-%!      0.42883;  
-%!      0.53010;  
-%!      0.80683 ];
-%!
 %! [Q,R] = qr(A);
-%! [Q,R] = qrinsert(Q,R,3,x);
+%! [Q,R] = qrinsert(Q,R,3,u);
 %! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - [A(:,1:2) x A(:,3)]),Inf) < norm(A)*1e1*eps)
-%! 
+%! assert(norm(vec(Q*R - [A(:,1:2) u A(:,3)]),Inf) < norm(A)*1e1*eps)
 %!test
-%! A = [0.620405 + 0.956953i  0.480013 + 0.048806i  0.402627 + 0.338171i;
-%!      0.589077 + 0.658457i  0.013205 + 0.279323i  0.229284 + 0.721929i;
-%!      0.092758 + 0.345687i  0.928679 + 0.241052i  0.764536 + 0.832406i;
-%!      0.912098 + 0.721024i  0.049018 + 0.269452i  0.730029 + 0.796517i;
-%!      0.112849 + 0.603871i  0.486352 + 0.142337i  0.355646 + 0.151496i ];
-%!
-%! x = [0.20351 + 0.05401i;
-%!      0.13141 + 0.43708i;
-%!      0.29808 + 0.08789i;
-%!      0.69821 + 0.38844i;
-%!      0.74871 + 0.25821i ];
-%!
-%! [Q,R] = qr(A);
-%! [Q,R] = qrinsert(Q,R,3,x);
+%! [Q,R] = qr(Ac);
+%! [Q,R] = qrinsert(Q,R,3,uc);
 %! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - [A(:,1:2) x A(:,3)]),Inf) < norm(A)*1e1*eps)
-%!
+%! assert(norm(vec(Q*R - [Ac(:,1:2) uc Ac(:,3)]),Inf) < norm(Ac)*1e1*eps)
 %!test
-%! A = [0.091364  0.613038  0.999083;
-%!      0.594638  0.425302  0.603537;
-%!      0.383594  0.291238  0.085574;
-%!      0.265712  0.268003  0.238409;
-%!      0.669966  0.743851  0.445057 ];
-%!
 %! x = [0.85082  0.76426  0.42883 ];
 %!
 %! [Q,R] = qr(A);
@@ -1020,21 +1086,43 @@
 %! assert(norm(vec(Q'*Q - eye(6)),Inf) < 1e1*eps)
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
 %! assert(norm(vec(Q*R - [A(1:2,:);x;A(3:5,:)]),Inf) < norm(A)*1e1*eps)
-%! 
 %!test
-%! A = [0.620405 + 0.956953i  0.480013 + 0.048806i  0.402627 + 0.338171i;
-%!      0.589077 + 0.658457i  0.013205 + 0.279323i  0.229284 + 0.721929i;
-%!      0.092758 + 0.345687i  0.928679 + 0.241052i  0.764536 + 0.832406i;
-%!      0.912098 + 0.721024i  0.049018 + 0.269452i  0.730029 + 0.796517i;
-%!      0.112849 + 0.603871i  0.486352 + 0.142337i  0.355646 + 0.151496i ];
-%!
 %! x = [0.20351 + 0.05401i  0.13141 + 0.43708i  0.29808 + 0.08789i ];
 %!
-%! [Q,R] = qr(A);
+%! [Q,R] = qr(Ac);
 %! [Q,R] = qrinsert(Q,R,3,x,'row');
 %! assert(norm(vec(Q'*Q - eye(6)),Inf) < 1e1*eps)
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - [A(1:2,:);x;A(3:5,:)]),Inf) < norm(A)*1e1*eps)
+%! assert(norm(vec(Q*R - [Ac(1:2,:);x;Ac(3:5,:)]),Inf) < norm(Ac)*1e1*eps)
+
+%!test
+%! [Q,R] = qr(single(A));
+%! [Q,R] = qrinsert(Q,R,3,single(u));
+%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - single([A(:,1:2) u A(:,3)])),Inf) < norm(single(A))*1e1*eps('single'))
+%!test
+%! [Q,R] = qr(single(Ac));
+%! [Q,R] = qrinsert(Q,R,3,single(uc));
+%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - single([Ac(:,1:2) uc Ac(:,3)])),Inf) < norm(single(Ac))*1e1*eps('single'))
+%!test
+%! x = single([0.85082  0.76426  0.42883 ]);
+%!
+%! [Q,R] = qr(single(A));
+%! [Q,R] = qrinsert(Q,R,3,x,'row');
+%! assert(norm(vec(Q'*Q - eye(6,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - single([A(1:2,:);x;A(3:5,:)])),Inf) < norm(single(A))*1e1*eps('single'))
+%!test
+%! x = single([0.20351 + 0.05401i  0.13141 + 0.43708i  0.29808 + 0.08789i ]);
+%!
+%! [Q,R] = qr(single(Ac));
+%! [Q,R] = qrinsert(Q,R,3,x,'row');
+%! assert(norm(vec(Q'*Q - eye(6,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - single([Ac(1:2,:);x;Ac(3:5,:)])),Inf) < norm(single(Ac))*1e1*eps('single'))
 */
 
 DEFUN_DLD (qrdelete, args, ,
@@ -1101,46 +1189,93 @@
 		    && argr.is_real_matrix ())
                   {
                     // real case
-                    Matrix Q = argq.matrix_value ();
-                    Matrix R = argr.matrix_value ();
+		    if (argq.is_single_type ()
+			|| argr.is_single_type ())
+		      {
+			FloatMatrix Q = argq.float_matrix_value ();
+			FloatMatrix R = argr.float_matrix_value ();
+
+			FloatQR fact (Q, R);
 
-                    QR fact (Q, R);
+			if (row) 
+			  fact.delete_row (j-1);
+			else 
+			  {
+			    fact.delete_col (j-1);
+
+			    if (! colp && k < m)
+			      fact.economize ();
+			  }
 
-                    if (row) 
-                      fact.delete_row (j-1);
-                    else 
-                      {
-                        fact.delete_col (j-1);
+			retval(1) = fact.R ();
+			retval(0) = fact.Q ();
+		      }
+		    else
+		      {
+			Matrix Q = argq.matrix_value ();
+			Matrix R = argr.matrix_value ();
+
+			QR fact (Q, R);
 
-                        if (! colp && k < m)
-                          fact.economize ();
-                      }
+			if (row) 
+			  fact.delete_row (j-1);
+			else 
+			  {
+			    fact.delete_col (j-1);
 
-                    retval(1) = fact.R ();
-                    retval(0) = fact.Q ();
+			    if (! colp && k < m)
+			      fact.economize ();
+			  }
+
+			retval(1) = fact.R ();
+			retval(0) = fact.Q ();
+		      }
                   }
                 else
                   {
                     // complex case
-                    ComplexMatrix Q = argq.complex_matrix_value ();
-                    ComplexMatrix R = argr.complex_matrix_value ();
+		    if (argq.is_single_type ()
+			|| argr.is_single_type ())
+		      {
+			FloatComplexMatrix Q = argq.float_complex_matrix_value ();
+			FloatComplexMatrix R = argr.float_complex_matrix_value ();
+
+			FloatComplexQR fact (Q, R);
 
-                    ComplexQR fact (Q, R);
+			if (row) 
+			  fact.delete_row (j-1);
+			else 
+			  {
+			    fact.delete_col (j-1);
+
+			    if (! colp && k < m)
+			      fact.economize ();
+			  }
 
-                    if (row) 
-                      fact.delete_row (j-1);
-                    else 
-                      {
-                        fact.delete_col (j-1);
+			retval(1) = fact.R ();
+			retval(0) = fact.Q ();
+		      }
+		    else
+		      {
+			ComplexMatrix Q = argq.complex_matrix_value ();
+			ComplexMatrix R = argr.complex_matrix_value ();
+
+			ComplexQR fact (Q, R);
 
-                        if (! colp && k < m)
-                          fact.economize ();
-                      }
+			if (row) 
+			  fact.delete_row (j-1);
+			else 
+			  {
+			    fact.delete_col (j-1);
 
-                    retval(1) = fact.R ();
-                    retval(0) = fact.Q ();
+			    if (! colp && k < m)
+			      fact.economize ();
+			  }
+
+			retval(1) = fact.R ();
+			retval(0) = fact.Q ();
+		      }
                   }
-
               }
             else
               error ("qrdelete: index j out of range");
@@ -1159,56 +1294,108 @@
  
 /*
 %!test
-%! A = [0.091364  0.613038  0.027504  0.999083;
-%!      0.594638  0.425302  0.562834  0.603537;
-%!      0.383594  0.291238  0.742073  0.085574;
-%!      0.265712  0.268003  0.783553  0.238409;
-%!      0.669966  0.743851  0.457255  0.445057 ];
+%! AA = [0.091364  0.613038  0.027504  0.999083;
+%!       0.594638  0.425302  0.562834  0.603537;
+%!       0.383594  0.291238  0.742073  0.085574;
+%!       0.265712  0.268003  0.783553  0.238409;
+%!       0.669966  0.743851  0.457255  0.445057 ];
 %!
-%! [Q,R] = qr(A);
+%! [Q,R] = qr(AA);
 %! [Q,R] = qrdelete(Q,R,3);
 %! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - [A(:,1:2) A(:,4)]),Inf) < norm(A)*1e1*eps)
+%! assert(norm(vec(Q*R - [AA(:,1:2) AA(:,4)]),Inf) < norm(AA)*1e1*eps)
 %! 
 %!test
-%! A = [0.364554 + 0.993117i  0.669818 + 0.510234i  0.426568 + 0.041337i  0.847051 + 0.233291i;
-%!      0.049600 + 0.242783i  0.448946 + 0.484022i  0.141155 + 0.074420i  0.446746 + 0.392706i;
-%!      0.581922 + 0.657416i  0.581460 + 0.030016i  0.219909 + 0.447288i  0.201144 + 0.069132i;
-%!      0.694986 + 0.000571i  0.682327 + 0.841712i  0.807537 + 0.166086i  0.192767 + 0.358098i;
-%!      0.945002 + 0.066788i  0.350492 + 0.642638i  0.579629 + 0.048102i  0.600170 + 0.636938i ] * I;
+%! AA = [0.364554 + 0.993117i  0.669818 + 0.510234i  0.426568 + 0.041337i  0.847051 + 0.233291i;
+%!       0.049600 + 0.242783i  0.448946 + 0.484022i  0.141155 + 0.074420i  0.446746 + 0.392706i;
+%!       0.581922 + 0.657416i  0.581460 + 0.030016i  0.219909 + 0.447288i  0.201144 + 0.069132i;
+%!       0.694986 + 0.000571i  0.682327 + 0.841712i  0.807537 + 0.166086i  0.192767 + 0.358098i;
+%!       0.945002 + 0.066788i  0.350492 + 0.642638i  0.579629 + 0.048102i  0.600170 + 0.636938i ] * I;
 %!
-%! [Q,R] = qr(A);
+%! [Q,R] = qr(AA);
 %! [Q,R] = qrdelete(Q,R,3);
 %! assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - [A(:,1:2) A(:,4)]),Inf) < norm(A)*1e1*eps)
+%! assert(norm(vec(Q*R - [AA(:,1:2) AA(:,4)]),Inf) < norm(AA)*1e1*eps)
 %!
 %!test
-%! A = [0.091364  0.613038  0.027504  0.999083;
-%!      0.594638  0.425302  0.562834  0.603537;
-%!      0.383594  0.291238  0.742073  0.085574;
-%!      0.265712  0.268003  0.783553  0.238409;
-%!      0.669966  0.743851  0.457255  0.445057 ];
+%! AA = [0.091364  0.613038  0.027504  0.999083;
+%!       0.594638  0.425302  0.562834  0.603537;
+%!       0.383594  0.291238  0.742073  0.085574;
+%!       0.265712  0.268003  0.783553  0.238409;
+%!       0.669966  0.743851  0.457255  0.445057 ];
 %!
-%! [Q,R] = qr(A);
+%! [Q,R] = qr(AA);
+%! [Q,R] = qrdelete(Q,R,3,'row');
+%! assert(norm(vec(Q'*Q - eye(4)),Inf) < 1e1*eps)
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - [AA(1:2,:);AA(4:5,:)]),Inf) < norm(AA)*1e1*eps)
+%! 
+%!test
+%! AA = [0.364554 + 0.993117i  0.669818 + 0.510234i  0.426568 + 0.041337i  0.847051 + 0.233291i;
+%!       0.049600 + 0.242783i  0.448946 + 0.484022i  0.141155 + 0.074420i  0.446746 + 0.392706i;
+%!       0.581922 + 0.657416i  0.581460 + 0.030016i  0.219909 + 0.447288i  0.201144 + 0.069132i;
+%!       0.694986 + 0.000571i  0.682327 + 0.841712i  0.807537 + 0.166086i  0.192767 + 0.358098i;
+%!       0.945002 + 0.066788i  0.350492 + 0.642638i  0.579629 + 0.048102i  0.600170 + 0.636938i ] * I;
+%!
+%! [Q,R] = qr(AA);
 %! [Q,R] = qrdelete(Q,R,3,'row');
 %! assert(norm(vec(Q'*Q - eye(4)),Inf) < 1e1*eps)
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - [A(1:2,:);A(4:5,:)]),Inf) < norm(A)*1e1*eps)
+%! assert(norm(vec(Q*R - [AA(1:2,:);AA(4:5,:)]),Inf) < norm(AA)*1e1*eps)
+
+%!test
+%! AA = single([0.091364  0.613038  0.027504  0.999083;
+%!              0.594638  0.425302  0.562834  0.603537;
+%!              0.383594  0.291238  0.742073  0.085574;
+%!              0.265712  0.268003  0.783553  0.238409;
+%!              0.669966  0.743851  0.457255  0.445057 ]);
+%!
+%! [Q,R] = qr(AA);
+%! [Q,R] = qrdelete(Q,R,3);
+%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - [AA(:,1:2) AA(:,4)]),Inf) < norm(AA)*1e1*eps('single'))
 %! 
 %!test
-%! A = [0.364554 + 0.993117i  0.669818 + 0.510234i  0.426568 + 0.041337i  0.847051 + 0.233291i;
-%!      0.049600 + 0.242783i  0.448946 + 0.484022i  0.141155 + 0.074420i  0.446746 + 0.392706i;
-%!      0.581922 + 0.657416i  0.581460 + 0.030016i  0.219909 + 0.447288i  0.201144 + 0.069132i;
-%!      0.694986 + 0.000571i  0.682327 + 0.841712i  0.807537 + 0.166086i  0.192767 + 0.358098i;
-%!      0.945002 + 0.066788i  0.350492 + 0.642638i  0.579629 + 0.048102i  0.600170 + 0.636938i ] * I;
+%! AA = single([0.364554 + 0.993117i  0.669818 + 0.510234i  0.426568 + 0.041337i  0.847051 + 0.233291i;
+%!              0.049600 + 0.242783i  0.448946 + 0.484022i  0.141155 + 0.074420i  0.446746 + 0.392706i;
+%!              0.581922 + 0.657416i  0.581460 + 0.030016i  0.219909 + 0.447288i  0.201144 + 0.069132i;
+%!              0.694986 + 0.000571i  0.682327 + 0.841712i  0.807537 + 0.166086i  0.192767 + 0.358098i;
+%!              0.945002 + 0.066788i  0.350492 + 0.642638i  0.579629 + 0.048102i  0.600170 + 0.636938i ]) * I;
+%!
+%! [Q,R] = qr(AA);
+%! [Q,R] = qrdelete(Q,R,3);
+%! assert(norm(vec(Q'*Q - eye(5,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - [AA(:,1:2) AA(:,4)]),Inf) < norm(AA)*1e1*eps('single'))
+%!
+%!test
+%! AA = single([0.091364  0.613038  0.027504  0.999083;
+%!              0.594638  0.425302  0.562834  0.603537;
+%!              0.383594  0.291238  0.742073  0.085574;
+%!              0.265712  0.268003  0.783553  0.238409;
+%!              0.669966  0.743851  0.457255  0.445057 ]);
 %!
-%! [Q,R] = qr(A);
+%! [Q,R] = qr(AA);
 %! [Q,R] = qrdelete(Q,R,3,'row');
-%! assert(norm(vec(Q'*Q - eye(4)),Inf) < 1e1*eps)
+%! assert(norm(vec(Q'*Q - eye(4,'single')),Inf) < 1e1*eps('single'))
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - [A(1:2,:);A(4:5,:)]),Inf) < norm(A)*1e1*eps)
+%! assert(norm(vec(Q*R - [AA(1:2,:);AA(4:5,:)]),Inf) < norm(AA)*1e1*eps('single'))
+%! 
+%!test
+%! AA = single([0.364554 + 0.993117i  0.669818 + 0.510234i  0.426568 + 0.041337i  0.847051 + 0.233291i;
+%!              0.049600 + 0.242783i  0.448946 + 0.484022i  0.141155 + 0.074420i  0.446746 + 0.392706i;
+%!              0.581922 + 0.657416i  0.581460 + 0.030016i  0.219909 + 0.447288i  0.201144 + 0.069132i;
+%!              0.694986 + 0.000571i  0.682327 + 0.841712i  0.807537 + 0.166086i  0.192767 + 0.358098i;
+%!              0.945002 + 0.066788i  0.350492 + 0.642638i  0.579629 + 0.048102i  0.600170 + 0.636938i ]) * I;
+%!
+%! [Q,R] = qr(AA);
+%! [Q,R] = qrdelete(Q,R,3,'row');
+%! assert(norm(vec(Q'*Q - eye(4,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - [AA(1:2,:);AA(4:5,:)]),Inf) < norm(AA)*1e1*eps('single'))
 */
 
 DEFUN_DLD (qrshift, args, ,
@@ -1256,26 +1443,56 @@
                   && argr.is_real_matrix ())
                 {
                   // all real case
-                  Matrix Q = argq.matrix_value ();
-                  Matrix R = argr.matrix_value ();
+		  if (argq.is_single_type () 
+		      && argr.is_single_type ())
+		    {
+		      FloatMatrix Q = argq.float_matrix_value ();
+		      FloatMatrix R = argr.float_matrix_value ();
+
+		      FloatQR fact (Q, R);
+		      fact.shift_cols (i-1, j-1);
 
-                  QR fact (Q, R);
-                  fact.shift_cols (i-1, j-1);
+		      retval(1) = fact.R ();
+		      retval(0) = fact.Q ();
+		    }
+		  else
+		    {
+		      Matrix Q = argq.matrix_value ();
+		      Matrix R = argr.matrix_value ();
 
-                  retval(1) = fact.R ();
-                  retval(0) = fact.Q ();
+		      QR fact (Q, R);
+		      fact.shift_cols (i-1, j-1);
+
+		      retval(1) = fact.R ();
+		      retval(0) = fact.Q ();
+		    }
                 }
               else
                 {
                   // complex case
-                  ComplexMatrix Q = argq.complex_matrix_value ();
-                  ComplexMatrix R = argr.complex_matrix_value ();
+		  if (argq.is_single_type () 
+		      && argr.is_single_type ())
+		    {
+		      FloatComplexMatrix Q = argq.float_complex_matrix_value ();
+		      FloatComplexMatrix R = argr.float_complex_matrix_value ();
 
-                  ComplexQR fact (Q, R);
-                  fact.shift_cols (i-1, j-1);
+		      FloatComplexQR fact (Q, R);
+		      fact.shift_cols (i-1, j-1);
                   
-                  retval(1) = fact.R ();
-                  retval(0) = fact.Q ();
+		      retval(1) = fact.R ();
+		      retval(0) = fact.Q ();
+		    }
+		  else
+		    {
+		      ComplexMatrix Q = argq.complex_matrix_value ();
+		      ComplexMatrix R = argr.complex_matrix_value ();
+
+		      ComplexQR fact (Q, R);
+		      fact.shift_cols (i-1, j-1);
+                  
+		      retval(1) = fact.R ();
+		      retval(0) = fact.Q ();
+		    }
                 }
             }
           else
@@ -1291,50 +1508,77 @@
 }
 /*
 %!test
-%! A = [0.091364  0.613038  0.999083;
-%!      0.594638  0.425302  0.603537;
-%!      0.383594  0.291238  0.085574;
-%!      0.265712  0.268003  0.238409;
-%!      0.669966  0.743851  0.445057 ].';
-%!
+%! AA = A.';
 %! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
 %!
-%! [Q,R] = qr(A);
+%! [Q,R] = qr(AA);
+%! [Q,R] = qrshift(Q,R,i,j);
+%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps)
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps)
+%! 
+%! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
+%!
+%! [Q,R] = qr(AA);
 %! [Q,R] = qrshift(Q,R,i,j);
 %! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps)
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - A(:,p)),Inf) < norm(A)*1e1*eps)
+%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps)
+%! 
+%!test
+%! AA = Ac.';
+%! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
+%!
+%! [Q,R] = qr(AA);
+%! [Q,R] = qrshift(Q,R,i,j);
+%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps)
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps)
 %! 
 %! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
 %!
-%! [Q,R] = qr(A);
+%! [Q,R] = qr(AA);
 %! [Q,R] = qrshift(Q,R,i,j);
 %! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps)
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - A(:,p)),Inf) < norm(A)*1e1*eps)
-%! 
+%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps)
+
+
 %!test
-%! A = [0.620405 + 0.956953i  0.480013 + 0.048806i  0.402627 + 0.338171i;
-%!      0.589077 + 0.658457i  0.013205 + 0.279323i  0.229284 + 0.721929i;
-%!      0.092758 + 0.345687i  0.928679 + 0.241052i  0.764536 + 0.832406i;
-%!      0.912098 + 0.721024i  0.049018 + 0.269452i  0.730029 + 0.796517i;
-%!      0.112849 + 0.603871i  0.486352 + 0.142337i  0.355646 + 0.151496i ].';
-%!
+%! AA = single (A).';
 %! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
 %!
-%! [Q,R] = qr(A);
+%! [Q,R] = qr(AA);
 %! [Q,R] = qrshift(Q,R,i,j);
-%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps)
+%! assert(norm(vec(Q'*Q - eye(3,'single')),Inf) < 1e1*eps('single'))
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - A(:,p)),Inf) < norm(A)*1e1*eps)
+%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single'))
 %! 
 %! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
 %!
-%! [Q,R] = qr(A);
+%! [Q,R] = qr(AA);
+%! [Q,R] = qrshift(Q,R,i,j);
+%! assert(norm(vec(Q'*Q - eye(3,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single'))
+%! 
+%!test
+%! AA = single(Ac).';
+%! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
+%!
+%! [Q,R] = qr(AA);
 %! [Q,R] = qrshift(Q,R,i,j);
-%! assert(norm(vec(Q'*Q - eye(3)),Inf) < 1e1*eps)
+%! assert(norm(vec(Q'*Q - eye(3,'single')),Inf) < 1e1*eps('single'))
 %! assert(norm(vec(triu(R)-R),Inf) == 0)
-%! assert(norm(vec(Q*R - A(:,p)),Inf) < norm(A)*1e1*eps)
+%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single'))
+%! 
+%! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
+%!
+%! [Q,R] = qr(AA);
+%! [Q,R] = qrshift(Q,R,i,j);
+%! assert(norm(vec(Q'*Q - eye(3,'single')),Inf) < 1e1*eps('single'))
+%! assert(norm(vec(triu(R)-R),Inf) == 0)
+%! assert(norm(vec(Q*R - AA(:,p)),Inf) < norm(AA)*1e1*eps('single'))
 */
 
 /*