changeset 2385:58386d13b8f1 octave-forge

Changed the directory structure of linear-algebra to match the package system
author hauberg
date Sun, 20 Aug 2006 14:58:16 +0000
parents 6131ff9d7bc2
children 19d651a0c71f
files main/linear-algebra/COPYING main/linear-algebra/CmplxGSVD.cc main/linear-algebra/CmplxGSVD.h main/linear-algebra/DESCRIPTION main/linear-algebra/GramSchmidt.cc main/linear-algebra/Makefile main/linear-algebra/dbleGSVD.cc main/linear-algebra/dbleGSVD.h main/linear-algebra/funm.m main/linear-algebra/gsvd.cc main/linear-algebra/inst/funm.m main/linear-algebra/inst/thfm.m main/linear-algebra/src/CmplxGSVD.cc main/linear-algebra/src/CmplxGSVD.h main/linear-algebra/src/GramSchmidt.cc main/linear-algebra/src/Makeconf.base main/linear-algebra/src/Makeconf.in main/linear-algebra/src/Makefile main/linear-algebra/src/autogen.sh main/linear-algebra/src/configure.base main/linear-algebra/src/dbleGSVD.cc main/linear-algebra/src/dbleGSVD.h main/linear-algebra/src/gsvd.cc main/linear-algebra/src/install-sh main/linear-algebra/src/octinst.sh.in main/linear-algebra/thfm.m
diffstat 26 files changed, 2777 insertions(+), 1346 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/COPYING	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,339 @@
+                    GNU GENERAL PUBLIC LICENSE
+                       Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+                          675 Mass Ave, Cambridge, MA 02139, USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The licenses for most software are designed to take away your
+freedom to share and change it.  By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users.  This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it.  (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.)  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+  To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have.  You must make sure that they, too, receive or can get the
+source code.  And you must show them these terms so they know their
+rights.
+
+  We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+  Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software.  If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+  Finally, any free program is threatened constantly by software
+patents.  We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary.  To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                    GNU GENERAL PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License.  The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language.  (Hereinafter, translation is included without limitation in
+the term "modification".)  Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope.  The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+  1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+  2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+    a) You must cause the modified files to carry prominent notices
+    stating that you changed the files and the date of any change.
+
+    b) You must cause any work that you distribute or publish, that in
+    whole or in part contains or is derived from the Program or any
+    part thereof, to be licensed as a whole at no charge to all third
+    parties under the terms of this License.
+
+    c) If the modified program normally reads commands interactively
+    when run, you must cause it, when started running for such
+    interactive use in the most ordinary way, to print or display an
+    announcement including an appropriate copyright notice and a
+    notice that there is no warranty (or else, saying that you provide
+    a warranty) and that users may redistribute the program under
+    these conditions, and telling the user how to view a copy of this
+    License.  (Exception: if the Program itself is interactive but
+    does not normally print such an announcement, your work based on
+    the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole.  If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works.  But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+  3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+    a) Accompany it with the complete corresponding machine-readable
+    source code, which must be distributed under the terms of Sections
+    1 and 2 above on a medium customarily used for software interchange; or,
+
+    b) Accompany it with a written offer, valid for at least three
+    years, to give any third party, for a charge no more than your
+    cost of physically performing source distribution, a complete
+    machine-readable copy of the corresponding source code, to be
+    distributed under the terms of Sections 1 and 2 above on a medium
+    customarily used for software interchange; or,
+
+    c) Accompany it with the information you received as to the offer
+    to distribute corresponding source code.  (This alternative is
+    allowed only for noncommercial distribution and only if you
+    received the program in object code or executable form with such
+    an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it.  For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable.  However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+  4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License.  Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+  5. You are not required to accept this License, since you have not
+signed it.  However, nothing else grants you permission to modify or
+distribute the Program or its derivative works.  These actions are
+prohibited by law if you do not accept this License.  Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+  6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions.  You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+  7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all.  For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices.  Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+  8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded.  In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+  9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number.  If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation.  If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+  10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission.  For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this.  Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+                            NO WARRANTY
+
+  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) 19yy  <name of author>
+
+    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; 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.
+
+    You should have received a copy of the GNU General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+    Gnomovision version 69, Copyright (C) 19yy name of author
+    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary.  Here is a sample; alter the names:
+
+  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+  `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+  <signature of Ty Coon>, 1 April 1989
+  Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs.  If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library.  If this is what you want to do, use the GNU Library General
+Public License instead of this License.
--- a/main/linear-algebra/CmplxGSVD.cc	Sun Aug 20 14:44:45 2006 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,262 +0,0 @@
-/*
-
-Copyright (C) 1996, 1997 John W. Eaton
-Copyright (C) 2006 Pascal Dupuis
-This file is part of Octave.
-
-Octave 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; either version 2, or (at your option) any
-later version.
-
-Octave 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.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, write to the Free
-Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
-02110-1301, USA.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <iostream>
-
-#include "CmplxGSVD.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (zggsvd, ZGGSVD)  
-   (
-     F77_CONST_CHAR_ARG_DECL, 	// JOBU    (input) CHARACTER*1
-     F77_CONST_CHAR_ARG_DECL, 	// JOBV    (input) CHARACTER*1
-     F77_CONST_CHAR_ARG_DECL, 	// JOBQ    (input) CHARACTER*1
-     const octave_idx_type&,	// M       (input) INTEGER
-     const octave_idx_type&,	// N       (input) INTEGER
-     const octave_idx_type&,	// P       (input) INTEGER
-     octave_idx_type &, 	// K       (output) INTEGER
-     octave_idx_type &,		// L       (output) INTEGER
-     Complex*,			// A       (input/output) COMPLEX*16 array, dimension (LDA,N)
-     const octave_idx_type&, 	// LDA     (input) INTEGER
-     Complex*, 			// B       (input/output) COMPLEX*16 array, dimension (LDB,N)
-     const octave_idx_type&, 	// LDB     (input) INTEGER
-     double*, 			// ALPHA   (output) DOUBLE PRECISION array, dimension (N)
-     double*, 			// BETA    (output) DOUBLE PRECISION array, dimension (N)
-     Complex*,			// U       (output) COMPLEX*16 array, dimension (LDU,M)
-     const octave_idx_type&,	// LDU     (input) INTEGER 
-     Complex*, 			// V       (output) COMPLEX*16 array, dimension (LDV,P)
-     const octave_idx_type&,	// LDV     (input) INTEGER
-     Complex*,			// Q       (output) COMPLEX*16 array, dimension (LDQ,N) 
-     const octave_idx_type&,	// LDQ     (input) INTEGER
-     Complex*, 			// WORK    (workspace) COMPLEX*16 array
-     double*,			// RWORK   (workspace) DOUBLE PRECISION array
-     int*,	 		// IWORK   (workspace/output) INTEGER array, dimension (N)
-     octave_idx_type&		// INFO    (output)INTEGER
-     F77_CHAR_ARG_LEN_DECL
-     F77_CHAR_ARG_LEN_DECL
-     F77_CHAR_ARG_LEN_DECL
-     );
-}
-
-ComplexMatrix
-ComplexGSVD::left_singular_matrix_A (void) const
-{
-  if (type_computed == GSVD::sigma_only)
-    {
-      (*current_liboctave_error_handler)
-	("dbleGSVD: U not computed because type == GSVD::sigma_only");
-      return ComplexMatrix ();
-    }
-  else
-    return left_smA;
-}
-
-ComplexMatrix
-ComplexGSVD::left_singular_matrix_B (void) const
-{
-  if (type_computed == GSVD::sigma_only)
-    {
-      (*current_liboctave_error_handler)
-	("dbleGSVD: V not computed because type == GSVD::sigma_only");
-      return ComplexMatrix ();
-    }
-  else
-    return left_smB;
-}
-
-ComplexMatrix
-ComplexGSVD::right_singular_matrix (void) const
-{
-  if (type_computed == GSVD::sigma_only)
-    {
-      (*current_liboctave_error_handler)
-	("dbleSVD: X not computed because type == GSVD::sigma_only");
-      return ComplexMatrix ();
-    }
-  else
-    return right_sm;
-}
-
-ComplexMatrix
-ComplexGSVD::R_matrix (void) const
-{
-  if (type_computed == GSVD::sigma_only)
-    {
-      (*current_liboctave_error_handler)
-	("dbleSVD: R not computed because type == GSVD::sigma_only");
-      return ComplexMatrix ();
-    }
-  else
-    return R;
-}
-
-octave_idx_type
-ComplexGSVD::init (const ComplexMatrix& a, const ComplexMatrix& b, 
-		   GSVD::type gsvd_type)
-{
-  octave_idx_type info;
-
-  octave_idx_type m = a.rows ();
-  octave_idx_type n = a.cols ();
-  octave_idx_type p = b.rows ();
-  
-  R = a;
-  Complex *tmp_dataA = R.fortran_vec ();
-  
-  ComplexMatrix btmp = b;
-  Complex *tmp_dataB = btmp.fortran_vec ();
-
-  // octave_idx_type min_mn = m < n ? m : n;
-
-  char jobu = 'U';
-  char jobv = 'V';
-  char jobq = 'Q';
-
-  octave_idx_type nrow_u = m;
-  octave_idx_type nrow_v = p;
-  octave_idx_type nrow_q = n;
-
-  octave_idx_type k, l;
-
-  switch (gsvd_type)
-    {
-
-    case GSVD::sigma_only:
-
-      // Note:  for this case, both jobu and jobv should be 'N', but
-      // there seems to be a bug in dgesvd from Lapack V2.0.  To
-      // demonstrate the bug, set both jobu and jobv to 'N' and find
-      // the singular values of [eye(3), eye(3)].  The result is
-      // [-sqrt(2), -sqrt(2), -sqrt(2)].
-      //
-      // For Lapack 3.0, this problem seems to be fixed.
-
-      jobu = 'N';
-      jobv = 'N';
-      jobq = 'N';
-      nrow_u = nrow_v = nrow_q = 1;
-      break;
-
-    default:
-      break;
-    }
-
-  type_computed = gsvd_type;
-
-  if (! (jobu == 'N' || jobu == 'O')) {
-    left_smA.resize (nrow_u, m);
-  }
-  
-  Complex *u = left_smA.fortran_vec ();
-
-  if (! (jobv == 'N' || jobv == 'O')) {
-    left_smB.resize (nrow_v, p);
-  }
-
-  Complex *v = left_smB.fortran_vec ();
-
-  sigmaA.resize (n, n);
-  // double *c_vec  = sigmaA.fortran_vec ();
-
-  sigmaB.resize (n, n);
-  // double *s_vec  = sigmaB.fortran_vec ();
-
-  if (! (jobq == 'N' || jobq == 'O')) {
-    right_sm.resize (nrow_q, n);
-  }
-  Complex *q = right_sm.fortran_vec ();  
-  
-  octave_idx_type lwork = 3*n;
-  lwork = lwork > m ? lwork : m;
-  lwork = (lwork > p ? lwork : p) + n;
-
-  Array<Complex>  work (lwork);
-  Array<double>   alpha (n);
-  Array<double>   beta (n);
-  Array<double>   rwork(2*n);
-  Array<int>      iwork (n);
-
-  F77_XFCN (zggsvd, ZGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
-			     F77_CONST_CHAR_ARG2 (&jobv, 1),
-			     F77_CONST_CHAR_ARG2 (&jobq, 1),
-			     m, n, p, k, l, tmp_dataA, m, 
-			     tmp_dataB, p, alpha.fortran_vec (),
-			     beta.fortran_vec (), u, nrow_u, 
-			     v, nrow_v, q, nrow_q, work.fortran_vec (), 
-			     rwork.fortran_vec (), iwork.fortran_vec (), 
-			     info
-			     F77_CHAR_ARG_LEN (1)
-			     F77_CHAR_ARG_LEN (1)
-			     F77_CHAR_ARG_LEN (1)));
-  
-  if (f77_exception_encountered)
-    (*current_liboctave_error_handler) ("unrecoverable error in zggsvd");
-
-  if (info < 0) {
-    (*current_liboctave_error_handler) ("zggsvd.f: argument %d illegal", -info);
-  }  else { 
-    if (info > 0) {
-      (*current_liboctave_error_handler) ("zggsvd.f: Jacobi-type procedure failed to converge."); 
-    } else {
-      octave_idx_type i;
-      for (i = 0; i < k; i++) {
-	sigmaA.xelem(i, i) = 1.0;
-	sigmaB.xelem(i, i) = 0.0;
-      }
-      if (m-k-l >= 0) { 
-	for (i = k; i < l; i++) {
-	  sigmaA.xelem(i, i) = alpha.elem(i);
-	  sigmaB.xelem(i, i) = beta.elem(i);
-	} 
-      } else {
-	for (i = k; i < m; i++) {
-	  sigmaA.xelem(i, i) = alpha.elem(i);
-	  sigmaB.xelem(i, i) = beta.elem(i);
-	}
-	for (i = k; i < m; i++) {
-          sigmaA.xelem(i, i) = alpha.elem(i);
-          sigmaB.xelem(i, i) = beta.elem(i);
-        }
-	for (i = m; i < k+l; i++) {
-          sigmaA.xelem(i, i) = 0.0;
-          sigmaB.xelem(i, i) = 1.0;;
-        }
-      }
-    }
-  }
-  return info;
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/main/linear-algebra/CmplxGSVD.h	Sun Aug 20 14:44:45 2006 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,106 +0,0 @@
-/*
-
-Copyright (C) 1996, 1997 John W. Eaton
-Copyright (C) 2006  Pascal Dupuis
-
-This file is part of Octave.
-
-Octave 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; either version 2, or (at your option) any
-later version.
-
-Octave 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.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, write to the Free
-Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
-02110-1301, USA.
-
-*/
-
-#if !defined (octave_ComplexGSVD_h)
-#define octave_ComplexGSVD_h 1
-
-#include <iostream>
-
-#include "dDiagMatrix.h"
-#include "CMatrix.h"
-#include "dbleGSVD.h"
-
-class
-ComplexGSVD
-{
-public:
-
-  ComplexGSVD (void) { }
-
-  ComplexGSVD (const ComplexMatrix& a, const ComplexMatrix& b,  
-	       GSVD::type gsvd_type = GSVD::economy) 
-    { 
-      init (a, b, gsvd_type); 
-    }
-
-  ComplexGSVD (const ComplexMatrix& a, const ComplexMatrix& b, 
-	       octave_idx_type& info, GSVD::type gsvd_type = GSVD::economy)
-    {
-      info = init (a, b, gsvd_type);
-    }
-
-  ComplexGSVD (const ComplexGSVD& a)
-    : type_computed (a.type_computed),
-      sigmaA (a.sigmaA), sigmaB (a.sigmaB), 
-      left_smA (a.left_smA), left_smB (a.left_smB), right_sm (a.right_sm),
-      R(a.R) { }
-
-  ComplexGSVD& operator = (const ComplexGSVD& a)
-    {
-      if (this != &a)
-	{
-	  type_computed = a.type_computed;
-	  sigmaA = a.sigmaA;
-	  sigmaB = a.sigmaB;
-	  left_smA = a.left_smA;
-	  left_smB = a.left_smB;
-	  right_sm = a.right_sm;
-	  R = a.R;
-	}
-
-      return *this;
-    }
-
-  ~ComplexGSVD (void) { }
-
-  DiagMatrix singular_values_A (void) const { return sigmaA; }
-  DiagMatrix singular_values_B (void) const { return sigmaB; }
-
-  ComplexMatrix left_singular_matrix_A (void) const;
-  ComplexMatrix left_singular_matrix_B (void) const;
-
-  ComplexMatrix right_singular_matrix (void) const;
-  ComplexMatrix R_matrix (void) const;
-
-  friend std::ostream&  operator << (std::ostream& os, const ComplexGSVD& a);
-
-private:
-
-  GSVD::type type_computed;
-
-  DiagMatrix sigmaA, sigmaB;
-  ComplexMatrix left_smA, left_smB;
-  ComplexMatrix right_sm, R;
-
-  octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b, 
-			GSVD::type gsvd_type = GSVD::economy);
-};
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/DESCRIPTION	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,11 @@
+Name: Linear-Algebra
+Version: 1.0.0
+Date: 2006-08-05
+Author: Various Authors
+Maintainer: The Octave Community
+Title: Linear Algebra.
+Description: Add a description to this package!
+Categories: Linear-Algebra
+Depends: octave (>= 2.9.7)
+License: GPL version 2 or later
+Url: http://octave.sf.net
--- a/main/linear-algebra/GramSchmidt.cc	Sun Aug 20 14:44:45 2006 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,87 +0,0 @@
-// Copyright (C) 2002 Andreas Stahel
-//
-// 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; 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.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
-
-#include <iostream>
-#include <cmath>
-
-#include <octave/oct.h>
-
-#include <octave/parse.h>
-#include <octave/pager.h>
-
-#define max(a,b) (((a)<(b)) ? (b) : (a))
-#define min(a,b) (((a)<(b)) ? (a) : (b))
-
-//////////////////////////////////////////////////
-
-
-void GramSchmidt(Matrix &V, ColumnVector &norms,int Vr, int Vc){
-  double tmp=0.0;
-  for(int i=0; i<Vc; i++){
-    tmp=0.0;  // normalize column i
-    for(int j=Vr-1;j>=0;j--){tmp += V(j,i)*V(j,i);}
-    tmp = norms(i)= sqrt(tmp);
-    for(int j=Vr-1;j>=0;j--){V(j,i)/=tmp;}
-    for(int k=i+1;k<Vc;k++){
-      tmp=0.0; // scalar product tmp=<V(:,i),V(:,k)>
-      for(int kk=Vr-1;kk>=0;kk--) tmp += V(kk,i)*V(kk,k);
-      // V(:,k) = V(:,k)-tmp*A(:,i)
-      for(int kk=Vr-1;kk>=0;kk--) V(kk,k) -= tmp*V(kk,i);  
-    };
-  };
-};
-
-//////////////////////////////////////////////////
-
-
-DEFUN_DLD (GramSchmidt, args, , "[...] = GramSchmidt(...)\n\
-  apply the Gram Schmidt reduction to the columns of a matrix V\n\
-\n\
-  Vout = GramSchmidt(V)\n\
-  [Vout, ColLength] = GramSchmidt(V)\n\
-\n\
-   V    is is a matrix of size mxn\n\
-   Vout is is a matrix of size mxn,\n\
-        the columns of Vout are orthonormalized and we have\n\
-        span(V(:,1:k)) = span(Vout(:,1:k)) for k=1...n\n\
-   ColLength is a vector containing the lengths of the column vectors of V\n\
-        during the Gram Schmidt algorithm\n\
-\n\
-   The implementation is based of the modified Gram Schmidt algorithm as\n\
-   described in \"Matrix Computations\" by G. Golub and C. van Loan")
-
-{
-  octave_value_list retval;
-  
-  int nargin = args.length ();
-  if (nargin != 1) {
-    print_usage ();
-    return retval;
-  }
-
-  octave_value V_arg = args(0);
-  int col = V_arg.columns();
-  int row = V_arg.rows();
-  Matrix V= V_arg.matrix_value();
-  
-  ColumnVector ColLength(col);
-  
-  GramSchmidt(V,ColLength,row,col);
-  
-  retval(0)= V;
-  retval(1)= ColLength;
-  return retval;
-}
--- a/main/linear-algebra/Makefile	Sun Aug 20 14:44:45 2006 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,48 +0,0 @@
-sinclude ../../Makeconf
-
-ifndef OCTAVE_FORGE
-MKOCTFILE = mkoctfile
-endif
-
-DEFINES = -DHAVE_CONFIG_H
-GSVD_OBJECTS = gsvd.o dbleGSVD.o CmplxGSVD.o 
-GSVD_TARGET = gsvd.oct
-OBJECTS = GramSchmidt.o $(GSVD_OBJECTS)
-TARGETS = $(GSVD_TARGET) GramSchmidt.oct
-
-MYDEPENDS = $(patsubst %.o,%.d,$(OBJECTS))
-
-ifeq ($(MAKECMDGOALS),all)
-  DEPENDS = $(MYDEPENDS)
-endif
-ifeq ($(MAKECMDGOALS),)
-  DEPENDS = $(MYDEPENDS)
-endif
-
-.PHONY: all clean count
-
-.SUFFIXES:
-
-.PRECIOUS: %.d %.o
-
-all : $(TARGETS)
-
-$(GSVD_TARGET) : $(DEPENDS) $(GSVD_OBJECTS)
-	$(MKOCTFILE) $(DEFINES) $(GSVD_OBJECTS) -o $@
-
-ifneq (,$(DEPENDS))
-  sinclude $(DEPENDS)
-endif
-
-%.d:%.cc
-	$(MKOCTFILE) $(DEFINES) -M $<
-
-%.o:%.cc
-%.o:%.cc %.d
-	$(MKOCTFILE) $(DEFINES) -c $<
-
-clean: 
-	rm -f $(TARGETS) $(MYDEPENDS) $(OBJECTS) *~ $(MYDEPENDS) octave-core
-
-count:
-	wc *{.cc,.h}
--- a/main/linear-algebra/dbleGSVD.cc	Sun Aug 20 14:44:45 2006 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,269 +0,0 @@
-/*
-
-Copyright (C) 1996, 1997 John W. Eaton
-Copyright (C) 2006 Pascal Dupuis
-This file is part of Octave.
-
-Octave 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; either version 2, or (at your option) any
-later version.
-
-Octave 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.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, write to the Free
-Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
-02110-1301, USA.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <iostream>
-
-#include "dbleGSVD.h"
-#include "f77-fcn.h"
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (dggsvd, DGGSVD)  
-   (
-     F77_CONST_CHAR_ARG_DECL, 	// JOBU    (input) CHARACTER*1
-     F77_CONST_CHAR_ARG_DECL, 	// JOBV    (input) CHARACTER*1
-     F77_CONST_CHAR_ARG_DECL, 	// JOBQ    (input) CHARACTER*1
-     const octave_idx_type&,	// M       (input) INTEGER
-     const octave_idx_type&,	// N       (input) INTEGER
-     const octave_idx_type&,	// P       (input) INTEGER
-     octave_idx_type &, 	// K       (output) INTEGER
-     octave_idx_type &,		// L       (output) INTEGER
-     double*,			// A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-     const octave_idx_type&, 	// LDA     (input) INTEGER
-     double*, 			// B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
-     const octave_idx_type&, 	// LDB     (input) INTEGER
-     double*, 			// ALPHA   (output) DOUBLE PRECISION array, dimension (N)
-     double*, 			// BETA    (output) DOUBLE PRECISION array, dimension (N)
-     double*,			// U       (output) DOUBLE PRECISION array, dimension (LDU,M)
-     const octave_idx_type&,	// LDU     (input) INTEGER 
-     double*, 			// V       (output) DOUBLE PRECISION array, dimension (LDV,P)
-     const octave_idx_type&,	// LDV     (input) INTEGER
-     double*,			// Q       (output) DOUBLE PRECISION array, dimension (LDQ,N) 
-     const octave_idx_type&,	// LDQ     (input) INTEGER
-     double*, 			// WORK    (workspace) DOUBLE PRECISION array
-     int*,	 		// IWORK   (workspace/output) INTEGER array, dimension (N)
-     octave_idx_type&		// INFO    (output)INTEGER
-     F77_CHAR_ARG_LEN_DECL
-     F77_CHAR_ARG_LEN_DECL
-     F77_CHAR_ARG_LEN_DECL
-     );
-}
-
-Matrix
-GSVD::left_singular_matrix_A (void) const
-{
-  if (type_computed == GSVD::sigma_only)
-    {
-      (*current_liboctave_error_handler)
-	("dbleGSVD: U not computed because type == GSVD::sigma_only");
-      return Matrix ();
-    }
-  else
-    return left_smA;
-}
-
-Matrix
-GSVD::left_singular_matrix_B (void) const
-{
-  if (type_computed == GSVD::sigma_only)
-    {
-      (*current_liboctave_error_handler)
-	("dbleGSVD: V not computed because type == GSVD::sigma_only");
-      return Matrix ();
-    }
-  else
-    return left_smB;
-}
-
-Matrix
-GSVD::right_singular_matrix (void) const
-{
-  if (type_computed == GSVD::sigma_only)
-    {
-      (*current_liboctave_error_handler)
-	("dbleSVD: X not computed because type == GSVD::sigma_only");
-      return Matrix ();
-    }
-  else
-    return right_sm;
-}
-Matrix
-GSVD::R_matrix (void) const
-{
-  if (type_computed == GSVD::sigma_only)
-    {
-      (*current_liboctave_error_handler)
-	("dbleSVD: R not computed because type == GSVD::sigma_only");
-      return Matrix ();
-    }
-  else
-    return R;
-}
-
-octave_idx_type
-GSVD::init (const Matrix& a, const Matrix& b, GSVD::type gsvd_type)
-{
-  octave_idx_type info;
-
-  octave_idx_type m = a.rows ();
-  octave_idx_type n = a.cols ();
-  octave_idx_type p = b.rows ();
-  
-  R = a;
-  double *tmp_dataA = R.fortran_vec ();
-  
-  Matrix btmp = b;
-  double *tmp_dataB = btmp.fortran_vec ();
-
-  // octave_idx_type min_mn = m < n ? m : n;
-
-  char jobu = 'U';
-  char jobv = 'V';
-  char jobq = 'Q';
-
-  octave_idx_type nrow_u = m;
-  octave_idx_type nrow_v = p;
-  octave_idx_type nrow_q = n;
-
-  octave_idx_type k, l;
-
-  switch (gsvd_type)
-    {
-
-    case GSVD::sigma_only:
-
-      // Note:  for this case, both jobu and jobv should be 'N', but
-      // there seems to be a bug in dgesvd from Lapack V2.0.  To
-      // demonstrate the bug, set both jobu and jobv to 'N' and find
-      // the singular values of [eye(3), eye(3)].  The result is
-      // [-sqrt(2), -sqrt(2), -sqrt(2)].
-      //
-      // For Lapack 3.0, this problem seems to be fixed.
-
-      jobu = 'N';
-      jobv = 'N';
-      jobq = 'N';
-      nrow_u = nrow_v = nrow_q = 1;
-      break;
-
-    default:
-      break;
-    }
-
-  type_computed = gsvd_type;
-
-  if (! (jobu == 'N' || jobu == 'O')) {
-    left_smA.resize (nrow_u, m);
-  }
-  
-  double *u = left_smA.fortran_vec ();
-
-  if (! (jobv == 'N' || jobv == 'O')) {
-    left_smB.resize (nrow_v, p);
-  }
-
-  double *v = left_smB.fortran_vec ();
-
-  sigmaA.resize (n, n);
-  // double *c_vec  = sigmaA.fortran_vec ();
-
-  sigmaB.resize (n, n);
-  // double *s_vec  = sigmaB.fortran_vec ();
-
-  if (! (jobq == 'N' || jobq == 'O')) {
-    right_sm.resize (nrow_q, n);
-  }
-  double *q = right_sm.fortran_vec ();  
-  
-  octave_idx_type lwork = 3*n;
-  lwork = lwork > m ? lwork : m;
-  lwork = (lwork > p ? lwork : p) + n;
-
-  Array<double> work (lwork);
-  Array<double> alpha (n);
-  Array<double> beta (n);
-  Array<int> 	iwork (n);
-
-  F77_XFCN (dggsvd, DGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
-			     F77_CONST_CHAR_ARG2 (&jobv, 1),
-			     F77_CONST_CHAR_ARG2 (&jobq, 1),
-			     m, n, p, k, l, tmp_dataA, m, 
-			     tmp_dataB, p, alpha.fortran_vec (),
-			     beta.fortran_vec (), u, nrow_u, 
-			     v, nrow_v, q, nrow_q, work.fortran_vec (), 
-			     iwork.fortran_vec (), info
-			     F77_CHAR_ARG_LEN (1)
-			     F77_CHAR_ARG_LEN (1)
-			     F77_CHAR_ARG_LEN (1)));
-  
-  if (f77_exception_encountered)
-    (*current_liboctave_error_handler) ("unrecoverable error in dggsvd");
- 
-  if (info < 0) {
-    (*current_liboctave_error_handler) ("dggsvd.f: argument %d illegal", -info);
-  } else {
-    if (info > 0) {
-      (*current_liboctave_error_handler) ("dggsvd.f: Jacobi-type procedure failed to converge."); 
-	
-    } else {
-      octave_idx_type i;
-      for (i = 0; i < k; i++) {
-	sigmaA.xelem(i, i) = 1.0;
-	sigmaB.xelem(i, i) = 0.0;
-      }
-      if (m-k-l >= 0) { 
-	for (i = k; i < l; i++) {
-	  sigmaA.xelem(i, i) = alpha.elem(i);
-	  sigmaB.xelem(i, i) = beta.elem(i);
-	} 
-      } else {
-	for (i = k; i < m; i++) {
-	  sigmaA.xelem(i, i) = alpha.elem(i);
-	  sigmaB.xelem(i, i) = beta.elem(i);
-	}
-	for (i = k; i < m; i++) {
-          sigmaA.xelem(i, i) = alpha.elem(i);
-          sigmaB.xelem(i, i) = beta.elem(i);
-        }
-	for (i = m; i < k+l; i++) {
-          sigmaA.xelem(i, i) = 0.0;
-          sigmaB.xelem(i, i) = 1.0;;
-        }
-      }
-    }
-  }
-  return info;
-}
-
-std::ostream&
-operator << (std::ostream& os, const GSVD& a)
-{
-  os << a.left_singular_matrix_A () << "\n";
-  os << a.left_singular_matrix_B () << "\n";
-  os << a.singular_values_A () << "\n";
-  os << a.singular_values_B () << "\n";
-  os << a.right_singular_matrix () << "\n";
-
-  return os;
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/main/linear-algebra/dbleGSVD.h	Sun Aug 20 14:44:45 2006 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,104 +0,0 @@
-/*
-
-Copyright (C) 1996, 1997 John W. Eaton
-Copyright (C) 2006  Pascal Dupuis
-
-This file is part of Octave.
-
-Octave 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; either version 2, or (at your option) any
-later version.
-
-Octave 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.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, write to the Free
-Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
-02110-1301, USA.
-
-*/
-
-#if !defined (octave_GSVD_h)
-#define octave_GSVD_h 1
-
-#include <iostream>
-
-#include "dDiagMatrix.h"
-#include "dMatrix.h"
-
-class
-GSVD
-{
-public:
-
-  enum type
-    {
-      economy,
-      sigma_only
-    };
-
-  GSVD (void) : sigmaA (), sigmaB (), left_smA (), left_smB (), right_sm () { }
-  GSVD (const Matrix& a, const Matrix& b, type gsvd_type = GSVD::economy) { init (a, b, gsvd_type); }
-
-  GSVD (const Matrix& a, const Matrix& b, octave_idx_type& info, type gsvd_type = GSVD::economy)
-    {
-      info = init (a, b, gsvd_type);
-    }
-
-  GSVD (const GSVD& a)
-    : type_computed (a.type_computed),
-      sigmaA (a.sigmaA), sigmaB (a.sigmaB), 
-      left_smA (a.left_smA), left_smB (a.left_smB), right_sm (a.right_sm),
-      R(a.R) { }
-
-  GSVD& operator = (const GSVD& a)
-    {
-      if (this != &a)
-	{
-	  type_computed = a.type_computed;
-	  sigmaA = a.sigmaA;
-	  sigmaB = a.sigmaB;
-	  left_smA = a.left_smA;
-	  left_smB = a.left_smB;
-	  right_sm = a.right_sm;
-	  R = a.R;
-	}
-
-      return *this;
-    }
-
-  ~GSVD (void) { }
-
-  DiagMatrix singular_values_A (void) const { return sigmaA; }
-  DiagMatrix singular_values_B (void) const { return sigmaB; }
-
-  Matrix left_singular_matrix_A (void) const;
-  Matrix left_singular_matrix_B (void) const;
-
-  Matrix right_singular_matrix (void) const;
-  Matrix R_matrix (void) const;
-
-  friend std::ostream&  operator << (std::ostream& os, const GSVD& a);
-
-private:
-
-  GSVD::type type_computed;
-
-  DiagMatrix sigmaA, sigmaB;
-  Matrix left_smA, left_smB;
-  Matrix right_sm, R;
-
-  octave_idx_type init (const Matrix& a, const Matrix& b, type gsvd_type = economy);
-};
-
-#endif
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- a/main/linear-algebra/funm.m	Sun Aug 20 14:44:45 2006 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,63 +0,0 @@
-## Copyright (C) 2000 P.R. Nienhuis
-##
-## 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; either version 2, 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.
-##
-## You should have received a copy of the GNU General Public License
-## along with this program; see the file COPYING.  If not, write to the
-## Free Software Foundation, 59 Temple Place - Suite 330, Boston, MA
-## 02111-1307, USA.
-
-## funm:  Matrix equivalent of function 'name'
-##
-## Usage:    B = funm(A, name)
-##  where    A = square non-singular matrix, provisionally
-##               real-valued
-##           B = square result matrix
-##        name = string, name of function to apply to A.
-##        args = any arguments to pass to function 'name'
-##               The function must accept a vector and apply
-##               element-wise to that vector.
-##
-## Example:  To compute sqrtm(A), you could use funm(A, 'sqrt')
-##
-## Note that you should not use funm for 'sqrt', 'log' or 'exp'; instead
-## use sqrtm, logm and expm which are more robust. Similarly,
-## trigonometric and hyperbolic functions (cos, sin, tan, cot, sec, csc,
-## cosh, sinh, tanh, coth, sech, csch) are better handled by thfm(A,
-## name), which defines them in terms of the more robust expm.
-
-## NOTE: the following comments are withheld until they can be verified
-##
-## If you have a network of coupled systems, where for the individual
-## systems a solution exists in terms of scalar variables, in many
-## cases the network might be solved using the same form of the
-## solution but with substituting the matrix equivalent of the function
-## applied to the scalar variables.
-## The approach is to do an eigen-analysis of the network to decouple
-## the systems, apply the scalar functions to the eigenvalues,
-## and then recombine the systems into a network.
-
-## Author: P.R. Nienhuis <106130.1515@compuserve.com>
-## Additions by P. Kienzle, .........
-## 2001-03-01 Paul Kienzle
-##    * generate error for repeated eigenvalues
-
-function B = funm(A, name)
-
-  if (nargin != 2 || !ischar(name) || ischar(A))
-    usage ("B = funm (A, 'f' [, args])");
-  endif
-
-  [V, D] = eig (A);
-  D = diag (feval (name, diag(D)));
-  B = V * D * inv (V);
-  
-endfunction
--- a/main/linear-algebra/gsvd.cc	Sun Aug 20 14:44:45 2006 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,269 +0,0 @@
-/*
-
-Copyright (C) 1996, 1997 John W. Eaton
-Copyright (C) 2006 Pascal Dupuis
-
-This file is part of Octave.
-
-Octave 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; either version 2, or (at your option) any
-later version.
-
-Octave 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.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, write to the Free
-Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
-02110-1301, USA.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "CmplxGSVD.h"
-#include "dbleGSVD.h"
-
-#include "defun-dld.h"
-#include "error.h"
-#include "gripes.h"
-#include "oct-obj.h"
-#include "pr-output.h"
-#include "utils.h"
-
-DEFUN_DLD (gsvd, args, nargout,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {@var{s} =} gsvd (@var{a}, @var{b})\n\
-@deftypefnx {Loadable Function} {[@var{u}, @var{v}, @var{c}, @var{s}, @var{x}] =} gsvd (@var{a}, @var{b})\n\
-@cindex generalised singular value decomposition\n\
-Compute the generalised singular value decomposition of (@var{a}, @var{b})\n\
-@iftex\n\
-@tex\n\
-$$\n\
- U^H A X = C R\n\
- V^H B X = S R\n\
- C*C + S*S = eye(columns(A))\n\
- R is upper triangular\n\
-$$\n\
-@end tex\n\
-@end iftex\n\
-@ifinfo\n\
-\n\
-@example\n\
-u' * a * x = c * r\n\
-v' * b * x = s * r'\n\
-c * c + s * s = eye(columns(a))\n\
-r is upper triangular\n\
-@end example\n\
-@end ifinfo\n\
-\n\
-The function @code{gsvd} normally returns the vector of singular values.\n\
-If asked for three return values, it computes\n\
-@iftex\n\
-@tex\n\
-$U$, $S$, and $V$.\n\
-@end tex\n\
-@end iftex\n\
-@ifinfo\n\
-U, S, and V.\n\
-@end ifinfo\n\
-For example,\n\
-\n\
-@example\n\
-svd (hilb (3))\n\
-@end example\n\
-\n\
-@noindent\n\
-returns\n\
-\n\
-@example\n\
-ans =\n\
-\n\
-  1.4083189\n\
-  0.1223271\n\
-  0.0026873\n\
-@end example\n\
-\n\
-@noindent\n\
-and\n\
-\n\
-@example\n\
-[u, s, v] = svd (hilb (3))\n\
-@end example\n\
-\n\
-@noindent\n\
-returns\n\
-\n\
-@example\n\
-u =\n\
-\n\
-  -0.82704   0.54745   0.12766\n\
-  -0.45986  -0.52829  -0.71375\n\
-  -0.32330  -0.64901   0.68867\n\
-\n\
-s =\n\
-\n\
-  1.40832  0.00000  0.00000\n\
-  0.00000  0.12233  0.00000\n\
-  0.00000  0.00000  0.00269\n\
-\n\
-v =\n\
-\n\
-  -0.82704   0.54745   0.12766\n\
-  -0.45986  -0.52829  -0.71375\n\
-  -0.32330  -0.64901   0.68867\n\
-@end example\n\
-\n\
-If given a second argument, @code{svd} returns an economy-sized\n\
-decomposition, eliminating the unnecessary rows or columns of @var{u} or\n\
-@var{v}.\n\
-@end deftypefn")
-{
-  octave_value_list retval;
-
-  int nargin = args.length ();
-
-  if (nargin < 2 || nargin > 2 || (nargout > 1 && (nargout < 5 || nargout > 6)))
-    {
-      print_usage ();
-      return retval;
-    }
-
-  octave_value argA = args(0), argB = args(1);
-
-  octave_idx_type nr = argA.rows ();
-  octave_idx_type nc = argA.columns ();
-
-  octave_idx_type  nn = argB.rows ();
-  octave_idx_type  np = argB.columns ();
-  
-  if (nr == 0 || nc == 0)
-    {
-      if (nargout >= 5)
-	{ 
-	  for (int i = 3; i <= nargout; i++)
-	    retval(i) = identity_matrix (nr, nr);
-	  retval(2) = Matrix (nr, nc);
-	  retval(1) = identity_matrix (nc, nc);
-	  retval(0) = identity_matrix (nc, nc);
-	}
-      else
-	retval(0) = Matrix (0, 1);
-    }
-  else
-    {
-      if ((nc != np) || (nn != np))
-	{
-	  print_usage ();
-	  return retval;
-	}
-
-      GSVD::type type = ((nargout == 0 || nargout == 1)
-			? GSVD::sigma_only
-			: GSVD::economy );
-
-      if (argA.is_real_type () && argB.is_real_type ())
-	{
-	  Matrix tmpA = argA.matrix_value ();
-	  Matrix tmpB = argB.matrix_value ();
-
-	  if (! error_state)
-	    {
-	      if (tmpA.any_element_is_inf_or_nan ())
-		{
-		  error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); 
-		  return retval;
-		}
-	      
-	      if (tmpB.any_element_is_inf_or_nan ())
-		{
-		  error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); 
-		  return retval;
-		}
-	      
-
-	      GSVD result (tmpA, tmpB, type);
-
-	      // DiagMatrix sigma = result.singular_values ();
-
-	      if (nargout == 0 || nargout == 1)
-		{
-		  DiagMatrix sigA =  result.singular_values_A ();
-		  DiagMatrix sigB =  result.singular_values_B ();
-		  for (int i = 0; i < nc; i++)
-		    tmpA.xelem(i, i) /= tmpB.xelem(i, i);
-		  retval(0) = sigA.diag();
-		}
-	      else
-		{ 
-		  if (nargout > 5) retval(5) = result.R_matrix ();
-		  retval(4) = result.right_singular_matrix ();
-		  retval(3) = result.singular_values_B ();
-		  retval(2) = result.singular_values_A ();
-		  retval(1) = result.left_singular_matrix_B ();
-		  retval(0) = result.left_singular_matrix_A ();
-		}
-	    }
-	}
-      else if (argA.is_complex_type () || argB.is_complex_type ())
-	{
-	  ComplexMatrix ctmpA = argA.complex_matrix_value ();
-	  ComplexMatrix ctmpB = argB.complex_matrix_value ();
-
-	  if (! error_state)
-	    {
-	      if (ctmpA.any_element_is_inf_or_nan ())
-		{
-		  error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); 
-		  return retval;
-		}
-	      if (ctmpB.any_element_is_inf_or_nan ())
-		{
-		  error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); 
-		  return retval;
-		}
-
-	      ComplexGSVD result (ctmpA, ctmpB, type);
-
-	      // DiagMatrix sigma = result.singular_values ();
-
-	      if (nargout == 0 || nargout == 1)
-		{
-		  DiagMatrix sigA =  result.singular_values_A ();
-		  DiagMatrix sigB =  result.singular_values_B ();
-		  for (int i = 0; i < nc; i++)
-		    ctmpA.xelem(i, i) /= ctmpB.xelem(i, i);
-		  retval(0) = sigA.diag();
-		}
-	      else
-		{
-		  if (nargout > 5) retval(5) = result.R_matrix ();
-		  retval(4) = result.right_singular_matrix ();
-		  retval(3) = result.singular_values_B ();
-		  retval(2) = result.singular_values_A ();
-		  retval(1) = result.left_singular_matrix_B ();
-		  retval(0) = result.left_singular_matrix_A ();
-		}
-	    }
-	}
-      else
-	{
-	  gripe_wrong_type_arg ("gsvd", argA);
-	  gripe_wrong_type_arg ("gsvd", argB);
-	  return retval;
-	}
-    }
-
-  return retval;
-}
-
-/*
-;;; Local Variables: ***
-;;; mode: C++ ***
-;;; End: ***
-*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/funm.m	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,63 @@
+## Copyright (C) 2000 P.R. Nienhuis
+##
+## 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; either version 2, 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.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; see the file COPYING.  If not, write to the
+## Free Software Foundation, 59 Temple Place - Suite 330, Boston, MA
+## 02111-1307, USA.
+
+## funm:  Matrix equivalent of function 'name'
+##
+## Usage:    B = funm(A, name)
+##  where    A = square non-singular matrix, provisionally
+##               real-valued
+##           B = square result matrix
+##        name = string, name of function to apply to A.
+##        args = any arguments to pass to function 'name'
+##               The function must accept a vector and apply
+##               element-wise to that vector.
+##
+## Example:  To compute sqrtm(A), you could use funm(A, 'sqrt')
+##
+## Note that you should not use funm for 'sqrt', 'log' or 'exp'; instead
+## use sqrtm, logm and expm which are more robust. Similarly,
+## trigonometric and hyperbolic functions (cos, sin, tan, cot, sec, csc,
+## cosh, sinh, tanh, coth, sech, csch) are better handled by thfm(A,
+## name), which defines them in terms of the more robust expm.
+
+## NOTE: the following comments are withheld until they can be verified
+##
+## If you have a network of coupled systems, where for the individual
+## systems a solution exists in terms of scalar variables, in many
+## cases the network might be solved using the same form of the
+## solution but with substituting the matrix equivalent of the function
+## applied to the scalar variables.
+## The approach is to do an eigen-analysis of the network to decouple
+## the systems, apply the scalar functions to the eigenvalues,
+## and then recombine the systems into a network.
+
+## Author: P.R. Nienhuis <106130.1515@compuserve.com>
+## Additions by P. Kienzle, .........
+## 2001-03-01 Paul Kienzle
+##    * generate error for repeated eigenvalues
+
+function B = funm(A, name)
+
+  if (nargin != 2 || !ischar(name) || ischar(A))
+    usage ("B = funm (A, 'f' [, args])");
+  endif
+
+  [V, D] = eig (A);
+  D = diag (feval (name, diag(D)));
+  B = V * D * inv (V);
+  
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/thfm.m	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,138 @@
+%USAGE  y = thfm ( x, MODE )
+%
+%       trigonometric/hyperbolic functions of square matrix x
+%
+%MODE	cos   sin   tan   sec   csc   cot
+%	cosh  sinh  tanh  sech  csch  coth
+%       acos  asin  atan  asec  acsc  acot
+%       acosh asinh atanh asech acsch acoth
+%       sqrt  log   exp
+%
+%NOTE		--- IMPORTANT ---
+%	This algorithm does  NOT USE an eigensystem
+%	similarity transformation. It maps the MODE
+%	functions to functions of expm, logm and sqrtm, 
+%       which are known to be robust with respect to
+%	non-diagonalizable ('defective') x
+%
+%EXA	thfm( x ,'cos' )  calculates  matrix cosinus
+%	EVEN IF input matrix x IS NOT DIAGONALIZABLE
+%
+%ASSOC	expm, sqrtm, logm, funm
+% Copyright	(C) 2001 Rolf Fabian <fabian@tu-cottbus.de>
+% 010213
+%	published under current GNU GENERAL PUBLIC LICENSE
+
+% 2001-03-15 Paul Kienzle
+%     * extend with inverse functions and power functions
+%     * optimize handling of real input
+
+function y=thfm(x,M)
+				#% minimal arg check only
+  if	nargin~=2||~ischar(M)||ischar(x)	
+    usage ("y = thfm (x, MODE)");
+  endif
+
+  ## look for known functions of sqrt, log, exp
+  I = eye(size(x));
+  match = 1;
+  len =  length(M);
+  if	len==3
+    
+    if	M=='cos',  
+      if (isreal(x))     y = real( expm( i*x ) );
+      else               y = ( expm( i*x ) + expm( -i*x ) ) / 2;
+      endif
+      
+    elseif	M=='sin',
+      if (isreal(x))     y = imag( expm( i*x ) );
+      else               y = ( expm( i*x ) - expm( -i*x ) ) / (2*i);
+      endif
+      
+    elseif	M=='tan',
+      if (isreal(x))     t = expm( i*x );    y = imag(t)/real(t);
+      else     	         t = expm( -2*i*x ); y = -i*(I-t)/(I+t);
+      endif
+      
+    elseif	M=='cot',		% == cos/sin
+      if (isreal(x))     t = expm( i*x );    y = real(t)/imag(t);
+      else	         t = expm( -2*i*x ); y = i*(I+t)/(I-t);
+      endif
+      
+    elseif	M=='sec',  
+      if (isreal(x))     y = inv( real(expm(i*x)) );
+      else               y = inv( expm(i*x)+expm(-i*x) )*2 ;
+      endif
+      
+    elseif	M=='csc',  
+      if (isreal(x))     y = inv( imag(expm(i*x)) );
+      else               y = inv( expm(i*x)-expm(-i*x) )*2i;
+      endif
+
+    elseif    M=='log',  y = logm(x);
+      
+    elseif    M=='exp',  y = expm(x);
+      
+    else match = 0;
+
+    endif
+    
+  elseif	len==4
+    
+    if      M=='cosh',   y = ( expm(x)+expm(-x) )/2;
+      
+    elseif  M=='sinh',   y = ( expm(x)-expm(-x) )/2;
+      
+    elseif  M=='tanh'    t = expm( -2*x ); y = (I - t)/(I + t);
+     
+    elseif  M=='coth', 	 t = expm( -2*x ); y = (I + t)/(I - t);
+      
+    elseif  M=='sech',   y = 2 * inv( expm(x) + expm(-x) );
+      
+    elseif  M=='csch',   y = 2 * inv( expm(x) - expm(-x) );
+      
+    elseif  M=='asin',   y = -i * logm( i*x + sqrtm(I - x*x) );
+      
+    elseif  M=='acos',   y =  i * logm( x - i*sqrtm(I - x*x) );
+
+    elseif  M=='atan',   y = -i/2 * logm( (I + i*x)/(I - i*x) );
+
+    elseif  M=='acot',   y =  i/2 * logm( (I + i*x)/(i*x - I) );
+
+    elseif  M=='asec',   y = i * logm( ( I - sqrtm(I - x*x) ) / x );
+
+    elseif  M=='acsc',   y = -i * logm( i*( I + sqrtm(I - x*x) ) / x );
+
+    elseif  M=='sqrt',   y = sqrtm(x);
+
+    else match = 0;
+
+    end
+
+  elseif   len==5
+
+    if      M=='asinh',  y = logm( x + sqrtm (x*x + I) );
+      
+    elseif  M=='acosh',  y = logm( x + sqrtm (x*x - I) );
+      
+    elseif  M=='atanh',  y = logm( (I + x)/(I - x) ) / 2;
+      
+    elseif  M=='acoth',  y = logm( (I + x)/(x - I) ) / 2;
+
+    elseif  M=='asech',  y = logm( (I + sqrtm (I - x*x)) / x );
+
+    elseif  M=='acsch',  y = logm( (I + sqrtm (I + x*x)) / x );
+
+    else match = 0;
+
+    endif
+
+  else match = 0;
+
+  endif
+
+  ## if no known function found, use generic solver
+  if (match == 0)
+    y = funm( x, M );
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/CmplxGSVD.cc	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,262 @@
+/*
+
+Copyright (C) 1996, 1997 John W. Eaton
+Copyright (C) 2006 Pascal Dupuis
+This file is part of Octave.
+
+Octave 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; either version 2, or (at your option) any
+later version.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <iostream>
+
+#include "CmplxGSVD.h"
+#include "f77-fcn.h"
+#include "lo-error.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (zggsvd, ZGGSVD)  
+   (
+     F77_CONST_CHAR_ARG_DECL, 	// JOBU    (input) CHARACTER*1
+     F77_CONST_CHAR_ARG_DECL, 	// JOBV    (input) CHARACTER*1
+     F77_CONST_CHAR_ARG_DECL, 	// JOBQ    (input) CHARACTER*1
+     const octave_idx_type&,	// M       (input) INTEGER
+     const octave_idx_type&,	// N       (input) INTEGER
+     const octave_idx_type&,	// P       (input) INTEGER
+     octave_idx_type &, 	// K       (output) INTEGER
+     octave_idx_type &,		// L       (output) INTEGER
+     Complex*,			// A       (input/output) COMPLEX*16 array, dimension (LDA,N)
+     const octave_idx_type&, 	// LDA     (input) INTEGER
+     Complex*, 			// B       (input/output) COMPLEX*16 array, dimension (LDB,N)
+     const octave_idx_type&, 	// LDB     (input) INTEGER
+     double*, 			// ALPHA   (output) DOUBLE PRECISION array, dimension (N)
+     double*, 			// BETA    (output) DOUBLE PRECISION array, dimension (N)
+     Complex*,			// U       (output) COMPLEX*16 array, dimension (LDU,M)
+     const octave_idx_type&,	// LDU     (input) INTEGER 
+     Complex*, 			// V       (output) COMPLEX*16 array, dimension (LDV,P)
+     const octave_idx_type&,	// LDV     (input) INTEGER
+     Complex*,			// Q       (output) COMPLEX*16 array, dimension (LDQ,N) 
+     const octave_idx_type&,	// LDQ     (input) INTEGER
+     Complex*, 			// WORK    (workspace) COMPLEX*16 array
+     double*,			// RWORK   (workspace) DOUBLE PRECISION array
+     int*,	 		// IWORK   (workspace/output) INTEGER array, dimension (N)
+     octave_idx_type&		// INFO    (output)INTEGER
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     );
+}
+
+ComplexMatrix
+ComplexGSVD::left_singular_matrix_A (void) const
+{
+  if (type_computed == GSVD::sigma_only)
+    {
+      (*current_liboctave_error_handler)
+	("dbleGSVD: U not computed because type == GSVD::sigma_only");
+      return ComplexMatrix ();
+    }
+  else
+    return left_smA;
+}
+
+ComplexMatrix
+ComplexGSVD::left_singular_matrix_B (void) const
+{
+  if (type_computed == GSVD::sigma_only)
+    {
+      (*current_liboctave_error_handler)
+	("dbleGSVD: V not computed because type == GSVD::sigma_only");
+      return ComplexMatrix ();
+    }
+  else
+    return left_smB;
+}
+
+ComplexMatrix
+ComplexGSVD::right_singular_matrix (void) const
+{
+  if (type_computed == GSVD::sigma_only)
+    {
+      (*current_liboctave_error_handler)
+	("dbleSVD: X not computed because type == GSVD::sigma_only");
+      return ComplexMatrix ();
+    }
+  else
+    return right_sm;
+}
+
+ComplexMatrix
+ComplexGSVD::R_matrix (void) const
+{
+  if (type_computed == GSVD::sigma_only)
+    {
+      (*current_liboctave_error_handler)
+	("dbleSVD: R not computed because type == GSVD::sigma_only");
+      return ComplexMatrix ();
+    }
+  else
+    return R;
+}
+
+octave_idx_type
+ComplexGSVD::init (const ComplexMatrix& a, const ComplexMatrix& b, 
+		   GSVD::type gsvd_type)
+{
+  octave_idx_type info;
+
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+  octave_idx_type p = b.rows ();
+  
+  R = a;
+  Complex *tmp_dataA = R.fortran_vec ();
+  
+  ComplexMatrix btmp = b;
+  Complex *tmp_dataB = btmp.fortran_vec ();
+
+  // octave_idx_type min_mn = m < n ? m : n;
+
+  char jobu = 'U';
+  char jobv = 'V';
+  char jobq = 'Q';
+
+  octave_idx_type nrow_u = m;
+  octave_idx_type nrow_v = p;
+  octave_idx_type nrow_q = n;
+
+  octave_idx_type k, l;
+
+  switch (gsvd_type)
+    {
+
+    case GSVD::sigma_only:
+
+      // Note:  for this case, both jobu and jobv should be 'N', but
+      // there seems to be a bug in dgesvd from Lapack V2.0.  To
+      // demonstrate the bug, set both jobu and jobv to 'N' and find
+      // the singular values of [eye(3), eye(3)].  The result is
+      // [-sqrt(2), -sqrt(2), -sqrt(2)].
+      //
+      // For Lapack 3.0, this problem seems to be fixed.
+
+      jobu = 'N';
+      jobv = 'N';
+      jobq = 'N';
+      nrow_u = nrow_v = nrow_q = 1;
+      break;
+
+    default:
+      break;
+    }
+
+  type_computed = gsvd_type;
+
+  if (! (jobu == 'N' || jobu == 'O')) {
+    left_smA.resize (nrow_u, m);
+  }
+  
+  Complex *u = left_smA.fortran_vec ();
+
+  if (! (jobv == 'N' || jobv == 'O')) {
+    left_smB.resize (nrow_v, p);
+  }
+
+  Complex *v = left_smB.fortran_vec ();
+
+  sigmaA.resize (n, n);
+  // double *c_vec  = sigmaA.fortran_vec ();
+
+  sigmaB.resize (n, n);
+  // double *s_vec  = sigmaB.fortran_vec ();
+
+  if (! (jobq == 'N' || jobq == 'O')) {
+    right_sm.resize (nrow_q, n);
+  }
+  Complex *q = right_sm.fortran_vec ();  
+  
+  octave_idx_type lwork = 3*n;
+  lwork = lwork > m ? lwork : m;
+  lwork = (lwork > p ? lwork : p) + n;
+
+  Array<Complex>  work (lwork);
+  Array<double>   alpha (n);
+  Array<double>   beta (n);
+  Array<double>   rwork(2*n);
+  Array<int>      iwork (n);
+
+  F77_XFCN (zggsvd, ZGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
+			     F77_CONST_CHAR_ARG2 (&jobv, 1),
+			     F77_CONST_CHAR_ARG2 (&jobq, 1),
+			     m, n, p, k, l, tmp_dataA, m, 
+			     tmp_dataB, p, alpha.fortran_vec (),
+			     beta.fortran_vec (), u, nrow_u, 
+			     v, nrow_v, q, nrow_q, work.fortran_vec (), 
+			     rwork.fortran_vec (), iwork.fortran_vec (), 
+			     info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+  
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in zggsvd");
+
+  if (info < 0) {
+    (*current_liboctave_error_handler) ("zggsvd.f: argument %d illegal", -info);
+  }  else { 
+    if (info > 0) {
+      (*current_liboctave_error_handler) ("zggsvd.f: Jacobi-type procedure failed to converge."); 
+    } else {
+      octave_idx_type i;
+      for (i = 0; i < k; i++) {
+	sigmaA.xelem(i, i) = 1.0;
+	sigmaB.xelem(i, i) = 0.0;
+      }
+      if (m-k-l >= 0) { 
+	for (i = k; i < l; i++) {
+	  sigmaA.xelem(i, i) = alpha.elem(i);
+	  sigmaB.xelem(i, i) = beta.elem(i);
+	} 
+      } else {
+	for (i = k; i < m; i++) {
+	  sigmaA.xelem(i, i) = alpha.elem(i);
+	  sigmaB.xelem(i, i) = beta.elem(i);
+	}
+	for (i = k; i < m; i++) {
+          sigmaA.xelem(i, i) = alpha.elem(i);
+          sigmaB.xelem(i, i) = beta.elem(i);
+        }
+	for (i = m; i < k+l; i++) {
+          sigmaA.xelem(i, i) = 0.0;
+          sigmaB.xelem(i, i) = 1.0;;
+        }
+      }
+    }
+  }
+  return info;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/CmplxGSVD.h	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,106 @@
+/*
+
+Copyright (C) 1996, 1997 John W. Eaton
+Copyright (C) 2006  Pascal Dupuis
+
+This file is part of Octave.
+
+Octave 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; either version 2, or (at your option) any
+later version.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+#if !defined (octave_ComplexGSVD_h)
+#define octave_ComplexGSVD_h 1
+
+#include <iostream>
+
+#include "dDiagMatrix.h"
+#include "CMatrix.h"
+#include "dbleGSVD.h"
+
+class
+ComplexGSVD
+{
+public:
+
+  ComplexGSVD (void) { }
+
+  ComplexGSVD (const ComplexMatrix& a, const ComplexMatrix& b,  
+	       GSVD::type gsvd_type = GSVD::economy) 
+    { 
+      init (a, b, gsvd_type); 
+    }
+
+  ComplexGSVD (const ComplexMatrix& a, const ComplexMatrix& b, 
+	       octave_idx_type& info, GSVD::type gsvd_type = GSVD::economy)
+    {
+      info = init (a, b, gsvd_type);
+    }
+
+  ComplexGSVD (const ComplexGSVD& a)
+    : type_computed (a.type_computed),
+      sigmaA (a.sigmaA), sigmaB (a.sigmaB), 
+      left_smA (a.left_smA), left_smB (a.left_smB), right_sm (a.right_sm),
+      R(a.R) { }
+
+  ComplexGSVD& operator = (const ComplexGSVD& a)
+    {
+      if (this != &a)
+	{
+	  type_computed = a.type_computed;
+	  sigmaA = a.sigmaA;
+	  sigmaB = a.sigmaB;
+	  left_smA = a.left_smA;
+	  left_smB = a.left_smB;
+	  right_sm = a.right_sm;
+	  R = a.R;
+	}
+
+      return *this;
+    }
+
+  ~ComplexGSVD (void) { }
+
+  DiagMatrix singular_values_A (void) const { return sigmaA; }
+  DiagMatrix singular_values_B (void) const { return sigmaB; }
+
+  ComplexMatrix left_singular_matrix_A (void) const;
+  ComplexMatrix left_singular_matrix_B (void) const;
+
+  ComplexMatrix right_singular_matrix (void) const;
+  ComplexMatrix R_matrix (void) const;
+
+  friend std::ostream&  operator << (std::ostream& os, const ComplexGSVD& a);
+
+private:
+
+  GSVD::type type_computed;
+
+  DiagMatrix sigmaA, sigmaB;
+  ComplexMatrix left_smA, left_smB;
+  ComplexMatrix right_sm, R;
+
+  octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b, 
+			GSVD::type gsvd_type = GSVD::economy);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/GramSchmidt.cc	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,87 @@
+// Copyright (C) 2002 Andreas Stahel
+//
+// 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; 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.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+#include <iostream>
+#include <cmath>
+
+#include <octave/oct.h>
+
+#include <octave/parse.h>
+#include <octave/pager.h>
+
+#define max(a,b) (((a)<(b)) ? (b) : (a))
+#define min(a,b) (((a)<(b)) ? (a) : (b))
+
+//////////////////////////////////////////////////
+
+
+void GramSchmidt(Matrix &V, ColumnVector &norms,int Vr, int Vc){
+  double tmp=0.0;
+  for(int i=0; i<Vc; i++){
+    tmp=0.0;  // normalize column i
+    for(int j=Vr-1;j>=0;j--){tmp += V(j,i)*V(j,i);}
+    tmp = norms(i)= sqrt(tmp);
+    for(int j=Vr-1;j>=0;j--){V(j,i)/=tmp;}
+    for(int k=i+1;k<Vc;k++){
+      tmp=0.0; // scalar product tmp=<V(:,i),V(:,k)>
+      for(int kk=Vr-1;kk>=0;kk--) tmp += V(kk,i)*V(kk,k);
+      // V(:,k) = V(:,k)-tmp*A(:,i)
+      for(int kk=Vr-1;kk>=0;kk--) V(kk,k) -= tmp*V(kk,i);  
+    };
+  };
+};
+
+//////////////////////////////////////////////////
+
+
+DEFUN_DLD (GramSchmidt, args, , "[...] = GramSchmidt(...)\n\
+  apply the Gram Schmidt reduction to the columns of a matrix V\n\
+\n\
+  Vout = GramSchmidt(V)\n\
+  [Vout, ColLength] = GramSchmidt(V)\n\
+\n\
+   V    is is a matrix of size mxn\n\
+   Vout is is a matrix of size mxn,\n\
+        the columns of Vout are orthonormalized and we have\n\
+        span(V(:,1:k)) = span(Vout(:,1:k)) for k=1...n\n\
+   ColLength is a vector containing the lengths of the column vectors of V\n\
+        during the Gram Schmidt algorithm\n\
+\n\
+   The implementation is based of the modified Gram Schmidt algorithm as\n\
+   described in \"Matrix Computations\" by G. Golub and C. van Loan")
+
+{
+  octave_value_list retval;
+  
+  int nargin = args.length ();
+  if (nargin != 1) {
+    print_usage ();
+    return retval;
+  }
+
+  octave_value V_arg = args(0);
+  int col = V_arg.columns();
+  int row = V_arg.rows();
+  Matrix V= V_arg.matrix_value();
+  
+  ColumnVector ColLength(col);
+  
+  GramSchmidt(V,ColLength,row,col);
+  
+  retval(0)= V;
+  retval(1)= ColLength;
+  return retval;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/Makeconf.base	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,92 @@
+
+## Makeconf is automatically generated from Makeconf.base and Makeconf.add
+## in the various subdirectories.  To regenerate, use ./autogen.sh to
+## create a new ./Makeconf.in, then use ./configure to generate a new
+## Makeconf.
+
+OCTAVE_FORGE = 1
+
+SHELL = @SHELL@
+
+canonical_host_type = @canonical_host_type@
+prefix = @prefix@
+exec_prefix = @exec_prefix@
+bindir = @bindir@
+mandir = @mandir@
+libdir = @libdir@
+datadir = @datadir@
+infodir = @infodir@
+includedir = @includedir@
+INSTALL = @INSTALL@
+INSTALL_PROGRAM = @INSTALL_PROGRAM@
+INSTALL_SCRIPT = @INSTALL_SCRIPT@
+INSTALL_DATA = @INSTALL_DATA@
+INSTALLOCT=octinst.sh
+
+DESTDIR =
+
+RANLIB = @RANLIB@
+STRIP = @STRIP@
+LN_S = @LN_S@
+MKOCTLINK = @MKOCTLINK@
+OCTLINK= @OCTLINK@
+
+AWK = @AWK@
+
+# Most octave programs will be compiled with $(MKOCTFILE).  Those which
+# cannot use mkoctfile directly can request the flags that mkoctfile 
+# would use as follows:
+#    FLAG = $(shell $(MKOCTFILE) -p FLAG)
+# The following flags are for compiling programs that are independent
+# of Octave.  How confusing.
+CC = @CC@
+CFLAGS = @CFLAGS@
+CPPFLAGS = @CPPFLAGS@
+CPICFLAG = @CPICFLAG@
+CXX = @CXX@
+CXXFLAGS = @CXXFLAGS@
+CXXPICFLAG = @CXXPICFLAG@
+F77 = @F77@
+FFLAGS = @FFLAGS@
+FPICFLAG = @FPICFLAG@
+
+OCTAVE = @OCTAVE@
+OCTAVE_VERSION = @OCTAVE_VERSION@
+MKOCTFILE = @MKOCTFILE@ -DHAVE_OCTAVE_$(ver) -v
+SHLEXT = @SHLEXT@
+
+@DEFHAVE_X@
+X_CFLAGS = @X_CFLAGS@
+X_LIBS = @X_LIBS@
+
+ver = @ver@
+MPATH = @mpath@
+OPATH = @opath@
+XPATH = @xpath@
+ALTMPATH = @altmpath@
+ALTOPATH = @altopath@
+
+HAVE_DO_FORTRAN_INDEXING = @HAVE_DO_FORTRAN_INDEXING@
+HAVE_PROPAGATE_EMPTY_MATRICES = @HAVE_PROPAGATE_EMPTY_MATRICES@
+HAVE_OK_TO_LOSE_IMAGINARY_PART = @HAVE_OK_TO_LOSE_IMAGINARY_PART@
+HAVE_ND_ARRAYS = @HAVE_ND_ARRAYS@
+TYPEID_HAS_CLASS = @TYPEID_HAS_CLASS@
+CLASS_HAS_LOAD_SAVE = @CLASS_HAS_LOAD_SAVE@
+HAVE_OCTAVE_MAP_INDEX = @HAVE_OCTAVE_MAP_INDEX@
+HAVE_OCTAVE_CONCAT = @HAVE_OCTAVE_CONCAT@
+HAVE_SWAP_BYTES = @HAVE_SWAP_BYTES@
+HAVE_OCTAVE_UPLUS = @HAVE_OCTAVE_UPLUS@
+
+MAKEINFO = @MAKEINFO@
+TEXI2DVI = @TEXI2DVI@
+TEXI2HTML = @TEXI2HTML@
+DVIPDF = @DVIPDF@
+DVIPS = @DVIPS@
+
+MKDOC = @MKDOC@
+MKTEXI = @MKTEXI@
+
+%.o: %.c ; $(MKOCTFILE) -c $<
+%.o: %.f ; $(MKOCTFILE) -c $<
+%.o: %.cc ; $(MKOCTFILE) -c $<
+%.oct: %.cc ; $(MKOCTFILE) $<
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/Makeconf.in	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,92 @@
+
+## Makeconf is automatically generated from Makeconf.base and Makeconf.add
+## in the various subdirectories.  To regenerate, use ./autogen.sh to
+## create a new ./Makeconf.in, then use ./configure to generate a new
+## Makeconf.
+
+OCTAVE_FORGE = 1
+
+SHELL = @SHELL@
+
+canonical_host_type = @canonical_host_type@
+prefix = @prefix@
+exec_prefix = @exec_prefix@
+bindir = @bindir@
+mandir = @mandir@
+libdir = @libdir@
+datadir = @datadir@
+infodir = @infodir@
+includedir = @includedir@
+INSTALL = @INSTALL@
+INSTALL_PROGRAM = @INSTALL_PROGRAM@
+INSTALL_SCRIPT = @INSTALL_SCRIPT@
+INSTALL_DATA = @INSTALL_DATA@
+INSTALLOCT=octinst.sh
+
+DESTDIR =
+
+RANLIB = @RANLIB@
+STRIP = @STRIP@
+LN_S = @LN_S@
+MKOCTLINK = @MKOCTLINK@
+OCTLINK= @OCTLINK@
+
+AWK = @AWK@
+
+# Most octave programs will be compiled with $(MKOCTFILE).  Those which
+# cannot use mkoctfile directly can request the flags that mkoctfile 
+# would use as follows:
+#    FLAG = $(shell $(MKOCTFILE) -p FLAG)
+# The following flags are for compiling programs that are independent
+# of Octave.  How confusing.
+CC = @CC@
+CFLAGS = @CFLAGS@
+CPPFLAGS = @CPPFLAGS@
+CPICFLAG = @CPICFLAG@
+CXX = @CXX@
+CXXFLAGS = @CXXFLAGS@
+CXXPICFLAG = @CXXPICFLAG@
+F77 = @F77@
+FFLAGS = @FFLAGS@
+FPICFLAG = @FPICFLAG@
+
+OCTAVE = @OCTAVE@
+OCTAVE_VERSION = @OCTAVE_VERSION@
+MKOCTFILE = @MKOCTFILE@ -DHAVE_OCTAVE_$(ver) -v
+SHLEXT = @SHLEXT@
+
+@DEFHAVE_X@
+X_CFLAGS = @X_CFLAGS@
+X_LIBS = @X_LIBS@
+
+ver = @ver@
+MPATH = @mpath@
+OPATH = @opath@
+XPATH = @xpath@
+ALTMPATH = @altmpath@
+ALTOPATH = @altopath@
+
+HAVE_DO_FORTRAN_INDEXING = @HAVE_DO_FORTRAN_INDEXING@
+HAVE_PROPAGATE_EMPTY_MATRICES = @HAVE_PROPAGATE_EMPTY_MATRICES@
+HAVE_OK_TO_LOSE_IMAGINARY_PART = @HAVE_OK_TO_LOSE_IMAGINARY_PART@
+HAVE_ND_ARRAYS = @HAVE_ND_ARRAYS@
+TYPEID_HAS_CLASS = @TYPEID_HAS_CLASS@
+CLASS_HAS_LOAD_SAVE = @CLASS_HAS_LOAD_SAVE@
+HAVE_OCTAVE_MAP_INDEX = @HAVE_OCTAVE_MAP_INDEX@
+HAVE_OCTAVE_CONCAT = @HAVE_OCTAVE_CONCAT@
+HAVE_SWAP_BYTES = @HAVE_SWAP_BYTES@
+HAVE_OCTAVE_UPLUS = @HAVE_OCTAVE_UPLUS@
+
+MAKEINFO = @MAKEINFO@
+TEXI2DVI = @TEXI2DVI@
+TEXI2HTML = @TEXI2HTML@
+DVIPDF = @DVIPDF@
+DVIPS = @DVIPS@
+
+MKDOC = @MKDOC@
+MKTEXI = @MKTEXI@
+
+%.o: %.c ; $(MKOCTFILE) -c $<
+%.o: %.f ; $(MKOCTFILE) -c $<
+%.o: %.cc ; $(MKOCTFILE) -c $<
+%.oct: %.cc ; $(MKOCTFILE) $<
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/Makefile	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,48 @@
+sinclude Makeconf
+
+ifndef OCTAVE_FORGE
+MKOCTFILE = mkoctfile
+endif
+
+DEFINES = -DHAVE_CONFIG_H
+GSVD_OBJECTS = gsvd.o dbleGSVD.o CmplxGSVD.o 
+GSVD_TARGET = gsvd.oct
+OBJECTS = GramSchmidt.o $(GSVD_OBJECTS)
+TARGETS = $(GSVD_TARGET) GramSchmidt.oct
+
+MYDEPENDS = $(patsubst %.o,%.d,$(OBJECTS))
+
+ifeq ($(MAKECMDGOALS),all)
+  DEPENDS = $(MYDEPENDS)
+endif
+ifeq ($(MAKECMDGOALS),)
+  DEPENDS = $(MYDEPENDS)
+endif
+
+.PHONY: all clean count
+
+.SUFFIXES:
+
+.PRECIOUS: %.d %.o
+
+all : $(TARGETS)
+
+$(GSVD_TARGET) : $(DEPENDS) $(GSVD_OBJECTS)
+	$(MKOCTFILE) $(DEFINES) $(GSVD_OBJECTS) -o $@
+
+ifneq (,$(DEPENDS))
+  sinclude $(DEPENDS)
+endif
+
+%.d:%.cc
+	$(MKOCTFILE) $(DEFINES) -M $<
+
+%.o:%.cc
+%.o:%.cc %.d
+	$(MKOCTFILE) $(DEFINES) -c $<
+
+clean: 
+	rm -f $(TARGETS) $(MYDEPENDS) $(OBJECTS) *~ $(MYDEPENDS) octave-core
+
+count:
+	wc *{.cc,.h}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/autogen.sh	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,33 @@
+#! /bin/sh
+
+## Generate ./configure
+rm -f configure.in
+echo "dnl --- DO NOT EDIT --- Automatically generated by autogen.sh" > configure.in
+cat configure.base >> configure.in
+files=`find . -name configure.add -print`
+if test ! -z "$files" ; then
+  cat $files >> configure.in
+fi
+cat <<EOF >> configure.in
+  AC_OUTPUT(\$CONFIGURE_OUTPUTS)
+  dnl XXX FIXME XXX chmod is not in autoconf's list of portable functions
+  chmod 0771 octinst.sh
+  echo " "
+  echo "  \"\\\$prefix\" is \$prefix"
+  echo "  \"\\\$exec_prefix\" is \$exec_prefix"
+  AC_MSG_RESULT([\$STATUS_MSG
+
+find . -name NOINSTALL -print    # shows which toolboxes won't be installed
+])
+EOF
+
+autoconf && rm -f configure.in
+
+## Generate ./Makeconf.in
+rm -f Makeconf.in
+cp Makeconf.base Makeconf.in
+files=`find . -name Makeconf.add -print`
+if test ! -z "$files" ; then
+  cat $files >> Makeconf.in
+fi
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/configure.base	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,520 @@
+dnl The configure script is generated by autogen.sh from configure.base 
+dnl and the various configure.add files in the source tree.  Edit 
+dnl configure.base and reprocess rather than modifying ./configure.
+
+dnl autoconf 2.13 certainly doesn't work! What is the minimum requirement?
+AC_PREREQ(2.2)
+
+AC_INIT(configure.base)
+
+PACKAGE=octave-forge
+MAJOR_VERSION=0
+MINOR_VERSION=1
+PATCH_LEVEL=0
+
+dnl Kill caching --- this ought to be the default
+define([AC_CACHE_LOAD], )dnl
+define([AC_CACHE_SAVE], )dnl
+
+dnl uncomment to put support files in another directory
+dnl AC_CONFIG_AUX_DIR(admin)
+
+VERSION=$MAJOR_VERSION.$MINOR_VERSION.$PATCH_LEVEL
+AC_SUBST(PACKAGE)
+AC_SUBST(VERSION)
+
+dnl need to find admin files, so keep track of the top dir.
+TOPDIR=`pwd`
+AC_SUBST(TOPDIR)
+
+dnl if mkoctfile doesn't work, then we need the following:
+dnl AC_PROG_CXX
+dnl AC_PROG_F77
+
+dnl Need C compiler regardless so define it in a way that
+dnl makes autoconf happy and we can override whatever we
+dnl need with mkoctfile -p.
+dnl XXX FIXME XXX should use mkoctfile to get CC and CFLAGS
+AC_PROG_CC
+
+dnl XXX FIXME XXX need tests for -p -c -s in mkoctfile.
+
+dnl *******************************************************************
+dnl Sort out mkoctfile version number and install paths
+
+dnl XXX FIXME XXX latest octave has octave-config so we don't
+dnl need to discover things here.  Doesn't have --exe-site-dir
+dnl but defines --oct-site-dir and --m-site-dir
+
+dnl Check for mkoctfile
+AC_CHECK_PROG(MKOCTFILE,mkoctfile,mkoctfile)
+test -z "$MKOCTFILE" &&	AC_MSG_WARN([no mkoctfile found on path])
+
+AC_SUBST(ver)
+AC_SUBST(subver)
+AC_SUBST(mpath)
+AC_SUBST(opath)
+AC_SUBST(xpath)
+AC_SUBST(altpath)
+AC_SUBST(altmpath)
+AC_SUBST(altopath)
+
+AC_ARG_WITH(path, 
+	[  --with-path             install path prefix],
+	[ path=$withval ])
+AC_ARG_WITH(mpath,
+	[  --with-mpath            override path for m-files],
+	[mpath=$withval])
+AC_ARG_WITH(opath,
+	[  --with-opath            override path for oct-files],
+	[opath=$withval])
+AC_ARG_WITH(xpath,
+	[  --with-xpath            override path for executables],
+	[xpath=$withval])
+AC_ARG_WITH(altpath, 
+	[  --with-altpath          alternative functions install path prefix],
+	[ altpath=$withval ])
+AC_ARG_WITH(altmpath,
+	[  --with-altmpath         override path for alternative m-files],
+	[altmpath=$withval])
+AC_ARG_WITH(altopath,
+	[  --with-altopath         override path for alternative oct-files],
+	[altopath=$withval])	
+
+if test -n "$path" ; then
+   test -z "$mpath" && mpath=$path 
+   test -z "$opath" && opath=$path/oct 
+   test -z "$xpath" && xpath=$path/bin
+   test -z "$altpath" && altpath=$path-alternatives
+fi
+
+if test -n "$altpath" ; then
+   test -z "$altmpath" && altmpath=$altpath 
+   test -z "$altopath" && altopath=$altpath/oct 
+fi
+
+dnl Don't query if path/ver are given in the configure environment
+#if test -z "$mpath" || test -z "$opath" || test -z "$xpath" || test -z "$altmpath" || test -z "$altopath" || test -z "$ver" ; then
+if test -z "$mpath" || test -z "$opath" || test -z "$xpath" || test -z "$ver" ; then
+   dnl Construct program to get mkoctfile version and local install paths
+   cat > conftest.cc <<EOF
+#include <octave/config.h>
+#include <octave/version.h>
+#include <octave/defaults.h>
+
+#define INFOV "\nINFOV=" OCTAVE_VERSION "\n"
+
+#define INFOH "\nINFOH=" OCTAVE_CANONICAL_HOST_TYPE "\n"
+
+#ifdef OCTAVE_LOCALVERFCNFILEDIR
+# define INFOM "\nINFOM=" OCTAVE_LOCALVERFCNFILEDIR "\n"
+#else
+# define INFOM "\nINFOM=" OCTAVE_LOCALFCNFILEPATH "\n"
+#endif
+
+#ifdef OCTAVE_LOCALVEROCTFILEDIR
+# define INFOO "\nINFOO=" OCTAVE_LOCALVEROCTFILEDIR "\n"
+#else
+# define INFOO "\nINFOO=" OCTAVE_LOCALOCTFILEPATH  "\n"
+#endif
+
+#ifdef OCTAVE_LOCALVERARCHLIBDIR
+# define INFOX "\nINFOX=" OCTAVE_LOCALVERARCHLIBDIR  "\n"
+#else
+# define INFOX "\nINFOX=" OCTAVE_LOCALARCHLIBDIR  "\n"
+#endif
+
+const char *infom = INFOM;
+const char *infoo = INFOO;
+const char *infox = INFOX;
+const char *infoh = INFOH;
+const char *infov = INFOV;
+EOF
+
+   dnl Compile program perhaps with a special version of mkoctfile
+   $MKOCTFILE conftest.cc || AC_MSG_ERROR(Could not run $MKOCTFILE)
+
+   dnl Strip the config info from the compiled file
+   eval `strings conftest.o | grep "^INFO.=" | sed -e "s,//.*$,,"`
+   rm -rf conftest*
+
+   dnl set the appropriate variables if they are not already set
+   ver=`echo $INFOV | sed -e "s/\.//" -e "s/\..*$//"`
+   subver=`echo $INFOV | sed -e "[s/^[^.]*[.][^.]*[.]//]"`
+   alt_mbase=`echo $INFOM | sed -e "[s,\/[^\/]*$,,]"`
+   alt_obase=`echo $INFOO | sed -e "[s,/site.*$,/site,]"`
+   test -z "$mpath" && mpath=$INFOM/octave-forge
+   test -z "$opath" && opath=$INFOO/octave-forge
+   test -z "$xpath" && xpath=$INFOX
+   test -z "$altmpath" && altmpath=$alt_mbase/octave-forge-alternatives/m
+   test -z "$altopath" && altopath=$alt_obase/octave-forge-alternatives/oct/$INFOH
+fi
+
+dnl *******************************************************************
+
+dnl XXX FIXME XXX Should we allow the user to override these?
+dnl Do we even need them?  The individual makefiles can call mkoctfile -p
+dnl themselves, so the only reason to keep them is for configure, and
+dnl for those things which are not built using mkoctfile (e.g., aurecord)
+dnl but it is not clear we should be using octave compile flags for those.
+
+dnl C compiler and flags
+AC_MSG_RESULT([retrieving compile and link flags from $MKOCTFILE])
+CC=`$MKOCTFILE -p CC`
+CFLAGS=`$MKOCTFILE -p CFLAGS`
+CPPFLAGS=`$MKOCTFILE -p CPPFLAGS`
+CPICFLAG=`$MKOCTFILE -p CPICFLAG`
+LDFLAGS=`$MKOCTFILE -p LDFLAGS`
+LIBS=`$MKOCTFILE -p LIBS`
+AC_SUBST(CC)
+AC_SUBST(CFLAGS)
+AC_SUBST(CPPFLAGS)
+AC_SUBST(CPICFLAG)
+
+dnl Fortran compiler and flags
+F77=`$MKOCTFILE -p F77`
+FFLAGS=`$MKOCTFILE -p FFLAGS`
+FPICFLAG=`$MKOCTFILE -p FPICFLAG`
+AC_SUBST(F77)
+AC_SUBST(FFLAGS)
+AC_SUBST(FPICFLAG)
+
+dnl C++ compiler and flags
+CXX=`$MKOCTFILE -p CXX`
+CXXFLAGS=`$MKOCTFILE -p CXXFLAGS`
+CXXPICFLAG=`$MKOCTFILE -p CXXPICFLAG`
+AC_SUBST(CXX)
+AC_SUBST(CXXFLAGS)
+AC_SUBST(CXXPICFLAG)
+
+dnl *******************************************************************
+
+dnl Check for features of your version of mkoctfile.
+dnl All checks should be designed so that the default
+dnl action if the tests are not performed is to do whatever
+dnl is appropriate for the most recent version of Octave.
+
+dnl Define the following macro:
+dnl    OF_CHECK_LIB(lib,fn,true,false,helpers)
+dnl This is just like AC_CHECK_LIB, but it doesn't update LIBS
+AC_DEFUN(OF_CHECK_LIB,
+[save_LIBS="$LIBS"
+AC_CHECK_LIB($1,$2,$3,$4,$5)
+LIBS="$save_LIBS"
+])
+
+dnl Define the following macro:
+dnl    TRY_MKOCTFILE(msg,program,action_if_true,action_if_false)
+dnl
+AC_DEFUN(TRY_MKOCTFILE,
+[AC_MSG_CHECKING($1)
+cat > conftest.cc << EOF
+#include <octave/config.h>
+$2
+EOF
+ac_try="$MKOCTFILE -c conftest.cc"
+if AC_TRY_EVAL(ac_try) ; then
+   AC_MSG_RESULT(yes)
+   $3
+else
+   AC_MSG_RESULT(no)
+   $4
+fi
+])
+
+dnl
+dnl Check if F77_FUNC works with MKOCTFILE
+dnl
+TRY_MKOCTFILE([for F77_FUNC],
+[int F77_FUNC (hello, HELLO) (const int &n);],,
+[MKOCTFILE="$MKOCTFILE -DF77_FUNC=F77_FCN"])
+
+dnl
+dnl Check if octave still uses SLList.h
+dnl
+TRY_MKOCTFILE([for SLList.h],[#include <octave/SLList.h>],
+[MKOCTFILE="$MKOCTFILE -DHAVE_SLLIST_H"],)
+
+dnl
+dnl Check if octave has lo_ieee_nan_value
+dnl
+TRY_MKOCTFILE([for lo_ieee_nan_value],
+[ #include <octave/lo-ieee.h>
+int test(void) { lo_ieee_nan_value(); }],,
+[MKOCTFILE="$MKOCTFILE -DUSE_OCTAVE_NAN"])
+
+dnl
+dnl Check if octave is needs octave_idx_type
+dnl
+TRY_MKOCTFILE([for octave_idx_type],
+[#include <octave/oct-types.h>
+octave_idx_type test(void) { octave_idx_type idx = 1; return idx; }],,
+[MKOCTFILE="$MKOCTFILE -Doctave_idx_type=int"])
+
+dnl
+dnl Check if octave uses quit.h
+dnl
+TRY_MKOCTFILE([for quit.h],[#include <octave/quit.h>],,
+[MKOCTFILE="$MKOCTFILE -DNEED_OCTAVE_QUIT"])
+
+dnl **********************************************************
+
+dnl Evaluate an expression in octave
+dnl
+dnl OCTAVE_EVAL(expr,var) -> var=expr
+dnl
+AC_DEFUN(OCTAVE_EVAL,
+[AC_MSG_CHECKING([for $1 in Octave])
+$2=`echo "disp($1)" | $OCTAVE -qf`
+AC_MSG_RESULT($$2)
+AC_SUBST($2)
+])
+
+dnl Check status of an octave variable
+dnl
+dnl OCTAVE_CHECK_EXIST(variable,action_if_true,action_if_false)
+dnl
+AC_DEFUN(OCTAVE_CHECK_EXIST,
+[AC_MSG_CHECKING([for $1 in Octave])
+if test `echo 'disp(exist("$1"))' | $OCTAVE -qf`X != 0X ; then
+   AC_MSG_RESULT(yes)
+   $2
+else
+   AC_MSG_RESULT(no)
+   $3
+fi
+])
+
+dnl should check that $(OCTAVE) --version matches $(MKOCTFILE) --version
+AC_CHECK_PROG(OCTAVE,octave,octave)
+OCTAVE_EVAL(OCTAVE_VERSION,OCTAVE_VERSION)
+
+dnl grab canonical host type so we can write system specific install stuff
+OCTAVE_EVAL(octave_config_info('canonical_host_type'),canonical_host_type)
+
+dnl grab SHLEXT from octave config
+OCTAVE_EVAL(octave_config_info('SHLEXT'),SHLEXT)
+
+AC_PROG_LN_S
+AC_PROG_INSTALL
+AC_PROG_RANLIB
+
+dnl Use $(COPY_FLAGS) to set options for cp when installing .oct files.
+COPY_FLAGS="-Rfp"
+case "$canonical_host_type" in
+    *-*-linux*)
+        COPY_FLAGS="-fdp"
+    ;;
+esac
+AC_SUBST(COPY_FLAGS)
+
+dnl Use $(STRIP) in the makefile to strip executables.  If not found, 
+dnl STRIP expands to ':', which in the makefile does nothing.
+dnl Don't need this for .oct files since mkoctfile handles them directly
+STRIP=${STRIP-strip}
+AC_CHECK_PROG(STRIP,$STRIP,$STRIP,:)
+
+dnl Strip on windows, don't strip on Mac OS/X or IRIX
+dnl For the rest, you can force strip using MKOCTFILE="mkoctfile -s"
+dnl or avoid strip using STRIP=: before ./configure
+case "$canonical_host_type" in
+    powerpc-apple-darwin*|*-sgi-*)
+	STRIP=:
+    ;;
+    *-cygwin-*|*-mingw-*) 
+	MKOCTFILE="$MKOCTFILE -s" 
+    ;;
+esac
+
+dnl Things needed to link to X11 programs
+dnl defines X_CFLAGS, X_LIBS
+AC_SUBST(DEFHAVE_X)
+AC_SUBST(X_LIBS)
+AC_SUBST(X_CFLAGS)
+AC_PATH_XTRA
+if test "$no_x" = yes ; then
+	DEFHAVE_X=
+	XSTATUS="no (plot/g{input,text,zoom,rab} will not work)"
+else
+	DEFHAVE_X="HAVE_X=1"
+	X_LIBS="$X_LIBS $X_PRE_LIBS -lX11 $X_EXTRA_LIBS"
+	XSTATUS="yes"
+fi
+
+OCTAVE_CHECK_EXIST(autoload,[
+	HAVE_AUTOLOAD="yes"
+	OCTLINK=.octlink
+	MKOCTLINK=$TOPDIR/admin/octlink.sh
+],[
+	HAVE_AUTOLOAD="no"
+	OCTLINK=.oct
+	MKOCTLINK=$LN_S
+])
+AC_SUBST(HAVE_AUTOLOAD)
+AC_SUBST(OCTLINK)
+AC_SUBST(MKOCTLINK)
+
+OCTAVE_CHECK_EXIST(do_fortran_indexing,
+	[HAVE_DO_FORTRAN_INDEXING="-DHAVE_DO_FORTRAN_INDEXING"],)
+AC_SUBST(HAVE_DO_FORTRAN_INDEXING)
+
+OCTAVE_CHECK_EXIST(propagate_empty_matrices,
+	[PROPAGATE_EMPTY_MATRICES="-DHAVE_PROPAGATE_EMPTY_MATRICES"],)
+AC_SUBST(HAVE_PROPAGATE_EMPTY_MATRICES)
+
+OCTAVE_CHECK_EXIST(ok_to_lose_imaginary_part,
+  [HAVE_OK_TO_LOSE_IMAGINARY_PART="-DHAVE_OK_TO_LOSE_IMAGINARY_PART"],)
+AC_SUBST(HAVE_OK_TO_LOSE_IMAGINARY_PART)
+
+dnl Test for N-dimensional Arrays
+TRY_MKOCTFILE([for N-dim arrays],
+[#include <octave/dim-vector.h>],
+[HAVE_ND_ARRAYS="-DHAVE_ND_ARRAYS"],)
+AC_SUBST(HAVE_ND_ARRAYS)
+
+OCTAVE_CHECK_EXIST(class,[TYPEID_HAS_CLASS="-DTYPEID_HAS_CLASS"],)
+AC_SUBST(TYPEID_HAS_CLASS)
+
+dnl Test for load/save functions in class
+TRY_MKOCTFILE([for load/save functions in class],
+[#include <octave/ov-scalar.h>
+int main (void) { octave_scalar a; a.load_ascii(std::cin); }],
+[CLASS_HAS_LOAD_SAVE="-DCLASS_HAS_LOAD_SAVE"],)
+AC_SUBST(CLASS_HAS_LOAD_SAVE)
+
+TRY_MKOCTFILE([for Octave_map indexing],
+[#include <octave/oct-map.h>
+int main(void) { Octave_map a; a[["key"]]; }],
+[HAVE_OCTAVE_MAP_INDEX="-DHAVE_OCTAVE_MAP_INDEX"],)
+AC_SUBST(HAVE_OCTAVE_MAP_INDEX)
+
+TRY_MKOCTFILE([for old Octave concatenation],
+[#include <octave/dNDArray.h>
+int main(void) { NDArray a(dim_vector(1,1)); Array<int> idx(2,0); a=concat(a,a,idx); }],
+[HAVE_OCTAVE_CONCAT="-DHAVE_OLD_OCTAVE_CONCAT"],)
+
+TRY_MKOCTFILE([for Octave concatenation],
+[#include <octave/dNDArray.h>
+int main(void) { NDArray a(dim_vector(1,1)); Array<int> idx(2,0); a=a.concat(a,idx); }],
+[HAVE_OCTAVE_CONCAT="-DHAVE_OCTAVE_CONCAT"],)
+AC_SUBST(HAVE_OCTAVE_CONCAT)
+
+TRY_MKOCTFILE([for swap_8_bytes],
+[#include <sys/types.h>
+#include <octave/config.h>
+#include <octave/byte-swap.h>
+int main(void) {long long a = 1; swap_8_bytes (&a,1);}],,
+[HAVE_SWAP_BYTES="-DHAVE_SWAP_BYTES"])
+AC_SUBST(HAVE_SWAP_BYTES)
+
+TRY_MKOCTFILE([for op_uplus],
+[#include <octave/config.h>
+#include <octave/ov.h>
+int main(void) {int i = octave_value::op_uplus;}],
+[HAVE_OCTAVE_UPLUS="-DHAVE_OCTAVE_UPLUS"],)
+AC_SUBST(HAVE_OCTAVE_UPLUS)
+
+dnl Test for the makeinfo program
+AC_CHECK_PROG(MAKEINFO,makeinfo,makeinfo)
+if [ test -n "$MAKEINFO" ]; then
+	dnl Check whether the makeinfo command accepts the 
+	dnl "--no-split" option
+	touch conftest.texi
+	AC_MSG_CHECKING([for makeinfo --no-split])
+	ac_try="$MAKEINFO --no-split conftest.texi"
+	if AC_TRY_EVAL(ac_try) ; then
+		MAKEINFO="$MAKEINFO --no-split"
+		AC_MSG_RESULT(yes)
+	else
+		AC_MSG_RESULT(no)
+	fi
+	rm -f conftest.*
+fi
+
+dnl Test for the texi2dvi program
+AC_CHECK_PROG(TEXI2DVI,texi2dvi,texi2dvi)
+if [ test -n "$TEXI2DVI" ]; then
+	dnl Check whether the texi2dvi command accepts the 
+	dnl "--clean" option
+	cat > conftest.texi <<EOF
+\input texinfo
+@bye
+EOF
+	AC_MSG_CHECKING([that texi2dvi runs])
+	ac_try="$TEXI2DVI conftest.texi > /dev/null"
+	if AC_TRY_EVAL(ac_try) ; then
+	    AC_MSG_RESULT(yes)
+	    AC_MSG_CHECKING([for texi2dvi --clean])
+	    ac_try="$TEXI2DVI --clean conftest.texi > /dev/null"
+	    if AC_TRY_EVAL(ac_try) ; then
+		TEXI2DVI="$TEXI2DVI --clean"
+		AC_MSG_RESULT(yes)
+	    else
+		AC_MSG_RESULT(no)
+	    fi
+	else
+	    TEXI2DVI=""
+	    AC_MSG_RESULT(no)
+	fi
+	rm -f conftest.*
+fi
+
+dnl Test for the texi2html program
+AC_CHECK_PROG(TEXI2HTML,texi2html,texi2html)
+if [ test -n "$TEXI2HTML" ]; then
+	STATUS="yes"
+	dnl Check whether the texi2html command accepts the 
+	dnl "-split_chapter -number" option
+	touch conftest.texi
+	AC_MSG_CHECKING([for texi2html --clean])
+	ac_try="$TEXI2HTML -split_chapter -number conftest.texi"
+	if AC_TRY_EVAL(ac_try) ; then
+		TEXI2HTML="$TEXI2HTML -split_chapter -number"
+		AC_MSG_RESULT(yes)
+	else
+		AC_MSG_RESULT(no)
+	fi
+	rm -f conftest.*
+        dnl TeTex 3.0 on Suse is leaving a conftest directory
+	rm -rf conftest
+fi
+
+dnl Test for the dvipdf program
+AC_CHECK_PROG(DVIPDF,dvipdf,dvipdf)
+
+dnl Test for the dvips program
+AC_CHECK_PROG(DVIPS,dvips,dvips)
+
+MKDOC=$TOPDIR/admin/mkdoc
+AC_SUBST(MKDOC)
+
+MKTEXI=$TOPDIR/admin/mktexi
+AC_SUBST(MKTEXI)
+
+CONFIGURE_OUTPUTS="Makeconf octinst.sh"
+STATUS_MSG="
+octave commands will install into the following directories:
+   m-files:   $mpath
+   oct-files: $opath
+   binaries:  $xpath
+alternatives:
+   m-files:   $altmpath
+   oct-files: $altopath
+
+shell commands will install into the following directories:
+   binaries:  $bindir
+   man pages: $mandir
+   libraries: $libdir
+   headers:   $includedir
+
+octave-forge is configured with
+   octave:      $OCTAVE (version $OCTAVE_VERSION)
+   mkoctfile:	$MKOCTFILE for Octave $subver
+   X11 support:	$XSTATUS
+   makeinfo:    $MAKEINFO
+   texi2dvi:    $TEXI2DVI
+   texi2html:   $TEXI2HTML
+   mkdoc:       $MKDOC
+   mktexi:      $MKTEXI
+   dvips:       $DVIPS
+   dvipdf:      $DVIPDF"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/dbleGSVD.cc	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,269 @@
+/*
+
+Copyright (C) 1996, 1997 John W. Eaton
+Copyright (C) 2006 Pascal Dupuis
+This file is part of Octave.
+
+Octave 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; either version 2, or (at your option) any
+later version.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <iostream>
+
+#include "dbleGSVD.h"
+#include "f77-fcn.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (dggsvd, DGGSVD)  
+   (
+     F77_CONST_CHAR_ARG_DECL, 	// JOBU    (input) CHARACTER*1
+     F77_CONST_CHAR_ARG_DECL, 	// JOBV    (input) CHARACTER*1
+     F77_CONST_CHAR_ARG_DECL, 	// JOBQ    (input) CHARACTER*1
+     const octave_idx_type&,	// M       (input) INTEGER
+     const octave_idx_type&,	// N       (input) INTEGER
+     const octave_idx_type&,	// P       (input) INTEGER
+     octave_idx_type &, 	// K       (output) INTEGER
+     octave_idx_type &,		// L       (output) INTEGER
+     double*,			// A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+     const octave_idx_type&, 	// LDA     (input) INTEGER
+     double*, 			// B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
+     const octave_idx_type&, 	// LDB     (input) INTEGER
+     double*, 			// ALPHA   (output) DOUBLE PRECISION array, dimension (N)
+     double*, 			// BETA    (output) DOUBLE PRECISION array, dimension (N)
+     double*,			// U       (output) DOUBLE PRECISION array, dimension (LDU,M)
+     const octave_idx_type&,	// LDU     (input) INTEGER 
+     double*, 			// V       (output) DOUBLE PRECISION array, dimension (LDV,P)
+     const octave_idx_type&,	// LDV     (input) INTEGER
+     double*,			// Q       (output) DOUBLE PRECISION array, dimension (LDQ,N) 
+     const octave_idx_type&,	// LDQ     (input) INTEGER
+     double*, 			// WORK    (workspace) DOUBLE PRECISION array
+     int*,	 		// IWORK   (workspace/output) INTEGER array, dimension (N)
+     octave_idx_type&		// INFO    (output)INTEGER
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     F77_CHAR_ARG_LEN_DECL
+     );
+}
+
+Matrix
+GSVD::left_singular_matrix_A (void) const
+{
+  if (type_computed == GSVD::sigma_only)
+    {
+      (*current_liboctave_error_handler)
+	("dbleGSVD: U not computed because type == GSVD::sigma_only");
+      return Matrix ();
+    }
+  else
+    return left_smA;
+}
+
+Matrix
+GSVD::left_singular_matrix_B (void) const
+{
+  if (type_computed == GSVD::sigma_only)
+    {
+      (*current_liboctave_error_handler)
+	("dbleGSVD: V not computed because type == GSVD::sigma_only");
+      return Matrix ();
+    }
+  else
+    return left_smB;
+}
+
+Matrix
+GSVD::right_singular_matrix (void) const
+{
+  if (type_computed == GSVD::sigma_only)
+    {
+      (*current_liboctave_error_handler)
+	("dbleSVD: X not computed because type == GSVD::sigma_only");
+      return Matrix ();
+    }
+  else
+    return right_sm;
+}
+Matrix
+GSVD::R_matrix (void) const
+{
+  if (type_computed == GSVD::sigma_only)
+    {
+      (*current_liboctave_error_handler)
+	("dbleSVD: R not computed because type == GSVD::sigma_only");
+      return Matrix ();
+    }
+  else
+    return R;
+}
+
+octave_idx_type
+GSVD::init (const Matrix& a, const Matrix& b, GSVD::type gsvd_type)
+{
+  octave_idx_type info;
+
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+  octave_idx_type p = b.rows ();
+  
+  R = a;
+  double *tmp_dataA = R.fortran_vec ();
+  
+  Matrix btmp = b;
+  double *tmp_dataB = btmp.fortran_vec ();
+
+  // octave_idx_type min_mn = m < n ? m : n;
+
+  char jobu = 'U';
+  char jobv = 'V';
+  char jobq = 'Q';
+
+  octave_idx_type nrow_u = m;
+  octave_idx_type nrow_v = p;
+  octave_idx_type nrow_q = n;
+
+  octave_idx_type k, l;
+
+  switch (gsvd_type)
+    {
+
+    case GSVD::sigma_only:
+
+      // Note:  for this case, both jobu and jobv should be 'N', but
+      // there seems to be a bug in dgesvd from Lapack V2.0.  To
+      // demonstrate the bug, set both jobu and jobv to 'N' and find
+      // the singular values of [eye(3), eye(3)].  The result is
+      // [-sqrt(2), -sqrt(2), -sqrt(2)].
+      //
+      // For Lapack 3.0, this problem seems to be fixed.
+
+      jobu = 'N';
+      jobv = 'N';
+      jobq = 'N';
+      nrow_u = nrow_v = nrow_q = 1;
+      break;
+
+    default:
+      break;
+    }
+
+  type_computed = gsvd_type;
+
+  if (! (jobu == 'N' || jobu == 'O')) {
+    left_smA.resize (nrow_u, m);
+  }
+  
+  double *u = left_smA.fortran_vec ();
+
+  if (! (jobv == 'N' || jobv == 'O')) {
+    left_smB.resize (nrow_v, p);
+  }
+
+  double *v = left_smB.fortran_vec ();
+
+  sigmaA.resize (n, n);
+  // double *c_vec  = sigmaA.fortran_vec ();
+
+  sigmaB.resize (n, n);
+  // double *s_vec  = sigmaB.fortran_vec ();
+
+  if (! (jobq == 'N' || jobq == 'O')) {
+    right_sm.resize (nrow_q, n);
+  }
+  double *q = right_sm.fortran_vec ();  
+  
+  octave_idx_type lwork = 3*n;
+  lwork = lwork > m ? lwork : m;
+  lwork = (lwork > p ? lwork : p) + n;
+
+  Array<double> work (lwork);
+  Array<double> alpha (n);
+  Array<double> beta (n);
+  Array<int> 	iwork (n);
+
+  F77_XFCN (dggsvd, DGGSVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
+			     F77_CONST_CHAR_ARG2 (&jobv, 1),
+			     F77_CONST_CHAR_ARG2 (&jobq, 1),
+			     m, n, p, k, l, tmp_dataA, m, 
+			     tmp_dataB, p, alpha.fortran_vec (),
+			     beta.fortran_vec (), u, nrow_u, 
+			     v, nrow_v, q, nrow_q, work.fortran_vec (), 
+			     iwork.fortran_vec (), info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+  
+  if (f77_exception_encountered)
+    (*current_liboctave_error_handler) ("unrecoverable error in dggsvd");
+ 
+  if (info < 0) {
+    (*current_liboctave_error_handler) ("dggsvd.f: argument %d illegal", -info);
+  } else {
+    if (info > 0) {
+      (*current_liboctave_error_handler) ("dggsvd.f: Jacobi-type procedure failed to converge."); 
+	
+    } else {
+      octave_idx_type i;
+      for (i = 0; i < k; i++) {
+	sigmaA.xelem(i, i) = 1.0;
+	sigmaB.xelem(i, i) = 0.0;
+      }
+      if (m-k-l >= 0) { 
+	for (i = k; i < l; i++) {
+	  sigmaA.xelem(i, i) = alpha.elem(i);
+	  sigmaB.xelem(i, i) = beta.elem(i);
+	} 
+      } else {
+	for (i = k; i < m; i++) {
+	  sigmaA.xelem(i, i) = alpha.elem(i);
+	  sigmaB.xelem(i, i) = beta.elem(i);
+	}
+	for (i = k; i < m; i++) {
+          sigmaA.xelem(i, i) = alpha.elem(i);
+          sigmaB.xelem(i, i) = beta.elem(i);
+        }
+	for (i = m; i < k+l; i++) {
+          sigmaA.xelem(i, i) = 0.0;
+          sigmaB.xelem(i, i) = 1.0;;
+        }
+      }
+    }
+  }
+  return info;
+}
+
+std::ostream&
+operator << (std::ostream& os, const GSVD& a)
+{
+  os << a.left_singular_matrix_A () << "\n";
+  os << a.left_singular_matrix_B () << "\n";
+  os << a.singular_values_A () << "\n";
+  os << a.singular_values_B () << "\n";
+  os << a.right_singular_matrix () << "\n";
+
+  return os;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/dbleGSVD.h	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,104 @@
+/*
+
+Copyright (C) 1996, 1997 John W. Eaton
+Copyright (C) 2006  Pascal Dupuis
+
+This file is part of Octave.
+
+Octave 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; either version 2, or (at your option) any
+later version.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+#if !defined (octave_GSVD_h)
+#define octave_GSVD_h 1
+
+#include <iostream>
+
+#include "dDiagMatrix.h"
+#include "dMatrix.h"
+
+class
+GSVD
+{
+public:
+
+  enum type
+    {
+      economy,
+      sigma_only
+    };
+
+  GSVD (void) : sigmaA (), sigmaB (), left_smA (), left_smB (), right_sm () { }
+  GSVD (const Matrix& a, const Matrix& b, type gsvd_type = GSVD::economy) { init (a, b, gsvd_type); }
+
+  GSVD (const Matrix& a, const Matrix& b, octave_idx_type& info, type gsvd_type = GSVD::economy)
+    {
+      info = init (a, b, gsvd_type);
+    }
+
+  GSVD (const GSVD& a)
+    : type_computed (a.type_computed),
+      sigmaA (a.sigmaA), sigmaB (a.sigmaB), 
+      left_smA (a.left_smA), left_smB (a.left_smB), right_sm (a.right_sm),
+      R(a.R) { }
+
+  GSVD& operator = (const GSVD& a)
+    {
+      if (this != &a)
+	{
+	  type_computed = a.type_computed;
+	  sigmaA = a.sigmaA;
+	  sigmaB = a.sigmaB;
+	  left_smA = a.left_smA;
+	  left_smB = a.left_smB;
+	  right_sm = a.right_sm;
+	  R = a.R;
+	}
+
+      return *this;
+    }
+
+  ~GSVD (void) { }
+
+  DiagMatrix singular_values_A (void) const { return sigmaA; }
+  DiagMatrix singular_values_B (void) const { return sigmaB; }
+
+  Matrix left_singular_matrix_A (void) const;
+  Matrix left_singular_matrix_B (void) const;
+
+  Matrix right_singular_matrix (void) const;
+  Matrix R_matrix (void) const;
+
+  friend std::ostream&  operator << (std::ostream& os, const GSVD& a);
+
+private:
+
+  GSVD::type type_computed;
+
+  DiagMatrix sigmaA, sigmaB;
+  Matrix left_smA, left_smB;
+  Matrix right_sm, R;
+
+  octave_idx_type init (const Matrix& a, const Matrix& b, type gsvd_type = economy);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/gsvd.cc	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,269 @@
+/*
+
+Copyright (C) 1996, 1997 John W. Eaton
+Copyright (C) 2006 Pascal Dupuis
+
+This file is part of Octave.
+
+Octave 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; either version 2, or (at your option) any
+later version.
+
+Octave 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.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CmplxGSVD.h"
+#include "dbleGSVD.h"
+
+#include "defun-dld.h"
+#include "error.h"
+#include "gripes.h"
+#include "oct-obj.h"
+#include "pr-output.h"
+#include "utils.h"
+
+DEFUN_DLD (gsvd, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{s} =} gsvd (@var{a}, @var{b})\n\
+@deftypefnx {Loadable Function} {[@var{u}, @var{v}, @var{c}, @var{s}, @var{x}] =} gsvd (@var{a}, @var{b})\n\
+@cindex generalised singular value decomposition\n\
+Compute the generalised singular value decomposition of (@var{a}, @var{b})\n\
+@iftex\n\
+@tex\n\
+$$\n\
+ U^H A X = C R\n\
+ V^H B X = S R\n\
+ C*C + S*S = eye(columns(A))\n\
+ R is upper triangular\n\
+$$\n\
+@end tex\n\
+@end iftex\n\
+@ifinfo\n\
+\n\
+@example\n\
+u' * a * x = c * r\n\
+v' * b * x = s * r'\n\
+c * c + s * s = eye(columns(a))\n\
+r is upper triangular\n\
+@end example\n\
+@end ifinfo\n\
+\n\
+The function @code{gsvd} normally returns the vector of singular values.\n\
+If asked for three return values, it computes\n\
+@iftex\n\
+@tex\n\
+$U$, $S$, and $V$.\n\
+@end tex\n\
+@end iftex\n\
+@ifinfo\n\
+U, S, and V.\n\
+@end ifinfo\n\
+For example,\n\
+\n\
+@example\n\
+svd (hilb (3))\n\
+@end example\n\
+\n\
+@noindent\n\
+returns\n\
+\n\
+@example\n\
+ans =\n\
+\n\
+  1.4083189\n\
+  0.1223271\n\
+  0.0026873\n\
+@end example\n\
+\n\
+@noindent\n\
+and\n\
+\n\
+@example\n\
+[u, s, v] = svd (hilb (3))\n\
+@end example\n\
+\n\
+@noindent\n\
+returns\n\
+\n\
+@example\n\
+u =\n\
+\n\
+  -0.82704   0.54745   0.12766\n\
+  -0.45986  -0.52829  -0.71375\n\
+  -0.32330  -0.64901   0.68867\n\
+\n\
+s =\n\
+\n\
+  1.40832  0.00000  0.00000\n\
+  0.00000  0.12233  0.00000\n\
+  0.00000  0.00000  0.00269\n\
+\n\
+v =\n\
+\n\
+  -0.82704   0.54745   0.12766\n\
+  -0.45986  -0.52829  -0.71375\n\
+  -0.32330  -0.64901   0.68867\n\
+@end example\n\
+\n\
+If given a second argument, @code{svd} returns an economy-sized\n\
+decomposition, eliminating the unnecessary rows or columns of @var{u} or\n\
+@var{v}.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+
+  int nargin = args.length ();
+
+  if (nargin < 2 || nargin > 2 || (nargout > 1 && (nargout < 5 || nargout > 6)))
+    {
+      print_usage ();
+      return retval;
+    }
+
+  octave_value argA = args(0), argB = args(1);
+
+  octave_idx_type nr = argA.rows ();
+  octave_idx_type nc = argA.columns ();
+
+  octave_idx_type  nn = argB.rows ();
+  octave_idx_type  np = argB.columns ();
+  
+  if (nr == 0 || nc == 0)
+    {
+      if (nargout >= 5)
+	{ 
+	  for (int i = 3; i <= nargout; i++)
+	    retval(i) = identity_matrix (nr, nr);
+	  retval(2) = Matrix (nr, nc);
+	  retval(1) = identity_matrix (nc, nc);
+	  retval(0) = identity_matrix (nc, nc);
+	}
+      else
+	retval(0) = Matrix (0, 1);
+    }
+  else
+    {
+      if ((nc != np) || (nn != np))
+	{
+	  print_usage ();
+	  return retval;
+	}
+
+      GSVD::type type = ((nargout == 0 || nargout == 1)
+			? GSVD::sigma_only
+			: GSVD::economy );
+
+      if (argA.is_real_type () && argB.is_real_type ())
+	{
+	  Matrix tmpA = argA.matrix_value ();
+	  Matrix tmpB = argB.matrix_value ();
+
+	  if (! error_state)
+	    {
+	      if (tmpA.any_element_is_inf_or_nan ())
+		{
+		  error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); 
+		  return retval;
+		}
+	      
+	      if (tmpB.any_element_is_inf_or_nan ())
+		{
+		  error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); 
+		  return retval;
+		}
+	      
+
+	      GSVD result (tmpA, tmpB, type);
+
+	      // DiagMatrix sigma = result.singular_values ();
+
+	      if (nargout == 0 || nargout == 1)
+		{
+		  DiagMatrix sigA =  result.singular_values_A ();
+		  DiagMatrix sigB =  result.singular_values_B ();
+		  for (int i = 0; i < nc; i++)
+		    tmpA.xelem(i, i) /= tmpB.xelem(i, i);
+		  retval(0) = sigA.diag();
+		}
+	      else
+		{ 
+		  if (nargout > 5) retval(5) = result.R_matrix ();
+		  retval(4) = result.right_singular_matrix ();
+		  retval(3) = result.singular_values_B ();
+		  retval(2) = result.singular_values_A ();
+		  retval(1) = result.left_singular_matrix_B ();
+		  retval(0) = result.left_singular_matrix_A ();
+		}
+	    }
+	}
+      else if (argA.is_complex_type () || argB.is_complex_type ())
+	{
+	  ComplexMatrix ctmpA = argA.complex_matrix_value ();
+	  ComplexMatrix ctmpB = argB.complex_matrix_value ();
+
+	  if (! error_state)
+	    {
+	      if (ctmpA.any_element_is_inf_or_nan ())
+		{
+		  error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); 
+		  return retval;
+		}
+	      if (ctmpB.any_element_is_inf_or_nan ())
+		{
+		  error ("gsvd: cannot take GSVD of matrix containing Inf or NaN values"); 
+		  return retval;
+		}
+
+	      ComplexGSVD result (ctmpA, ctmpB, type);
+
+	      // DiagMatrix sigma = result.singular_values ();
+
+	      if (nargout == 0 || nargout == 1)
+		{
+		  DiagMatrix sigA =  result.singular_values_A ();
+		  DiagMatrix sigB =  result.singular_values_B ();
+		  for (int i = 0; i < nc; i++)
+		    ctmpA.xelem(i, i) /= ctmpB.xelem(i, i);
+		  retval(0) = sigA.diag();
+		}
+	      else
+		{
+		  if (nargout > 5) retval(5) = result.R_matrix ();
+		  retval(4) = result.right_singular_matrix ();
+		  retval(3) = result.singular_values_B ();
+		  retval(2) = result.singular_values_A ();
+		  retval(1) = result.left_singular_matrix_B ();
+		  retval(0) = result.left_singular_matrix_A ();
+		}
+	    }
+	}
+      else
+	{
+	  gripe_wrong_type_arg ("gsvd", argA);
+	  gripe_wrong_type_arg ("gsvd", argB);
+	  return retval;
+	}
+    }
+
+  return retval;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/install-sh	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,251 @@
+#!/bin/sh
+#
+# install - install a program, script, or datafile
+# This comes from X11R5 (mit/util/scripts/install.sh).
+#
+# Copyright 1991 by the Massachusetts Institute of Technology
+#
+# Permission to use, copy, modify, distribute, and sell this software and its
+# documentation for any purpose is hereby granted without fee, provided that
+# the above copyright notice appear in all copies and that both that
+# copyright notice and this permission notice appear in supporting
+# documentation, and that the name of M.I.T. not be used in advertising or
+# publicity pertaining to distribution of the software without specific,
+# written prior permission.  M.I.T. makes no representations about the
+# suitability of this software for any purpose.  It is provided "as is"
+# without express or implied warranty.
+#
+# Calling this script install-sh is preferred over install.sh, to prevent
+# `make' implicit rules from creating a file called install from it
+# when there is no Makefile.
+#
+# This script is compatible with the BSD install script, but was written
+# from scratch.  It can only install one file at a time, a restriction
+# shared with many OS's install programs.
+
+
+# set DOITPROG to echo to test this script
+
+# Don't use :- since 4.3BSD and earlier shells don't like it.
+doit="${DOITPROG-}"
+
+
+# put in absolute paths if you don't have them in your path; or use env. vars.
+
+mvprog="${MVPROG-mv}"
+cpprog="${CPPROG-cp}"
+chmodprog="${CHMODPROG-chmod}"
+chownprog="${CHOWNPROG-chown}"
+chgrpprog="${CHGRPPROG-chgrp}"
+stripprog="${STRIPPROG-strip}"
+rmprog="${RMPROG-rm}"
+mkdirprog="${MKDIRPROG-mkdir}"
+
+transformbasename=""
+transform_arg=""
+instcmd="$mvprog"
+chmodcmd="$chmodprog 0755"
+chowncmd=""
+chgrpcmd=""
+stripcmd=""
+rmcmd="$rmprog -f"
+mvcmd="$mvprog"
+src=""
+dst=""
+dir_arg=""
+
+while [ x"$1" != x ]; do
+    case $1 in
+	-c) instcmd="$cpprog"
+	    shift
+	    continue;;
+
+	-d) dir_arg=true
+	    shift
+	    continue;;
+
+	-m) chmodcmd="$chmodprog $2"
+	    shift
+	    shift
+	    continue;;
+
+	-o) chowncmd="$chownprog $2"
+	    shift
+	    shift
+	    continue;;
+
+	-g) chgrpcmd="$chgrpprog $2"
+	    shift
+	    shift
+	    continue;;
+
+	-s) stripcmd="$stripprog"
+	    shift
+	    continue;;
+
+	-t=*) transformarg=`echo $1 | sed 's/-t=//'`
+	    shift
+	    continue;;
+
+	-b=*) transformbasename=`echo $1 | sed 's/-b=//'`
+	    shift
+	    continue;;
+
+	*)  if [ x"$src" = x ]
+	    then
+		src=$1
+	    else
+		# this colon is to work around a 386BSD /bin/sh bug
+		:
+		dst=$1
+	    fi
+	    shift
+	    continue;;
+    esac
+done
+
+if [ x"$src" = x ]
+then
+	echo "install:	no input file specified"
+	exit 1
+else
+	true
+fi
+
+if [ x"$dir_arg" != x ]; then
+	dst=$src
+	src=""
+	
+	if [ -d $dst ]; then
+		instcmd=:
+		chmodcmd=""
+	else
+		instcmd=mkdir
+	fi
+else
+
+# Waiting for this to be detected by the "$instcmd $src $dsttmp" command
+# might cause directories to be created, which would be especially bad 
+# if $src (and thus $dsttmp) contains '*'.
+
+	if [ -f $src -o -d $src ]
+	then
+		true
+	else
+		echo "install:  $src does not exist"
+		exit 1
+	fi
+	
+	if [ x"$dst" = x ]
+	then
+		echo "install:	no destination specified"
+		exit 1
+	else
+		true
+	fi
+
+# If destination is a directory, append the input filename; if your system
+# does not like double slashes in filenames, you may need to add some logic
+
+	if [ -d $dst ]
+	then
+		dst="$dst"/`basename $src`
+	else
+		true
+	fi
+fi
+
+## this sed command emulates the dirname command
+dstdir=`echo $dst | sed -e 's,[^/]*$,,;s,/$,,;s,^$,.,'`
+
+# Make sure that the destination directory exists.
+#  this part is taken from Noah Friedman's mkinstalldirs script
+
+# Skip lots of stat calls in the usual case.
+if [ ! -d "$dstdir" ]; then
+defaultIFS='	
+'
+IFS="${IFS-${defaultIFS}}"
+
+oIFS="${IFS}"
+# Some sh's can't handle IFS=/ for some reason.
+IFS='%'
+set - `echo ${dstdir} | sed -e 's@/@%@g' -e 's@^%@/@'`
+IFS="${oIFS}"
+
+pathcomp=''
+
+while [ $# -ne 0 ] ; do
+	pathcomp="${pathcomp}${1}"
+	shift
+
+	if [ ! -d "${pathcomp}" ] ;
+        then
+		$mkdirprog "${pathcomp}"
+	else
+		true
+	fi
+
+	pathcomp="${pathcomp}/"
+done
+fi
+
+if [ x"$dir_arg" != x ]
+then
+	$doit $instcmd $dst &&
+
+	if [ x"$chowncmd" != x ]; then $doit $chowncmd $dst; else true ; fi &&
+	if [ x"$chgrpcmd" != x ]; then $doit $chgrpcmd $dst; else true ; fi &&
+	if [ x"$stripcmd" != x ]; then $doit $stripcmd $dst; else true ; fi &&
+	if [ x"$chmodcmd" != x ]; then $doit $chmodcmd $dst; else true ; fi
+else
+
+# If we're going to rename the final executable, determine the name now.
+
+	if [ x"$transformarg" = x ] 
+	then
+		dstfile=`basename $dst`
+	else
+		dstfile=`basename $dst $transformbasename | 
+			sed $transformarg`$transformbasename
+	fi
+
+# don't allow the sed command to completely eliminate the filename
+
+	if [ x"$dstfile" = x ] 
+	then
+		dstfile=`basename $dst`
+	else
+		true
+	fi
+
+# Make a temp file name in the proper directory.
+
+	dsttmp=$dstdir/#inst.$$#
+
+# Move or copy the file name to the temp name
+
+	$doit $instcmd $src $dsttmp &&
+
+	trap "rm -f ${dsttmp}" 0 &&
+
+# and set any options; do chmod last to preserve setuid bits
+
+# If any of these fail, we abort the whole thing.  If we want to
+# ignore errors from any of these, just make sure not to ignore
+# errors from the above "$doit $instcmd $src $dsttmp" command.
+
+	if [ x"$chowncmd" != x ]; then $doit $chowncmd $dsttmp; else true;fi &&
+	if [ x"$chgrpcmd" != x ]; then $doit $chgrpcmd $dsttmp; else true;fi &&
+	if [ x"$stripcmd" != x ]; then $doit $stripcmd $dsttmp; else true;fi &&
+	if [ x"$chmodcmd" != x ]; then $doit $chmodcmd $dsttmp; else true;fi &&
+
+# Now rename the file to the real destination.
+
+	$doit $rmcmd -f $dstdir/$dstfile &&
+	$doit $mvcmd $dsttmp $dstdir/$dstfile 
+
+fi &&
+
+
+exit 0
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/src/octinst.sh.in	Sun Aug 20 14:58:16 2006 +0000
@@ -0,0 +1,93 @@
+#! /bin/sh
+
+# octinst.sh source mpath opath xpath [altmpath altopath]
+
+# Copies all m-files and oct-files from the source directory to the
+# mpath and opath respectively.  Preserves links.  Files in
+# source/data are copied to mpath.  Files in the source/bin are copied
+# to xpath. m-files and oct-files in source/alternatives are copied to 
+# altmpath and altopath respectively
+
+if test $# -lt 4 ; then
+    echo 'Not enough arguments'
+    exit 1
+fi
+
+# interpret input parameters
+source=$1; shift
+mpath=$1; shift
+opath=$1; shift
+xpath=$1; shift
+if test $# -ge 1; then altmpath=$1; shift; fi
+if test $# -ge 1; then altopath=$1; shift; fi
+INSTALL="@INSTALL@"
+INSTALL_DATA="@INSTALL_DATA@"
+INSTALL_PROGRAM="@INSTALL_PROGRAM@"
+INSTALL_SCRIPT="@INSTALL_SCRIPT@"
+MKPKGADD="@TOPDIR@/admin/mkpkgadd"
+COPY_FLAGS="@COPY_FLAGS@"
+
+# grab the m-files
+files=`echo $source/*.m`
+if test "$files" != "$source/*.m" ; then
+    $INSTALL -d $mpath
+    $INSTALL_DATA $files $mpath
+fi
+
+# grab the oct-files
+files=`echo $source/*.oct`
+if test "$files" != "$source/*.oct" ; then
+    $INSTALL -d $opath
+## Grrr... install doesn't preserve links.  Hope this works.
+    cp $COPY_FLAGS $files $opath
+fi
+
+# install alternatives
+if test -d "$source/alternatives" ; then
+    # m-files
+    files=`echo $source/alternatives/*.m`
+    if test "$files" != "$source/alternatives/*.m" ; then
+        $INSTALL -d $altmpath
+	$INSTALL_DATA $files $altmpath
+    fi
+    # oct-files
+    files=`echo $source/alternatives/*.oct`
+    if test "$files" != "$source/alternatives/*.oct" ; then
+        $INSTALL -d $altopath
+	$INSTALL_DATA $files $altopath
+    fi
+fi
+
+# Create PKG_ADD, and destroy it immediately if it is empty
+# XXX FIXME XXX no PKG_ADD created if only oct-files and no m-files.
+if test -d "$mpath" ; then
+    $MKPKGADD $source > $mpath/PKG_ADD
+    if test -z "`cat $mpath/PKG_ADD`" ; then rm -f $mpath/PKG_ADD;  fi
+fi
+# PKG_ADD for alternatives
+if test -d "$source/alternatives" -a -d "$altmpath" ; then
+    $MKPKGADD $source/alternatives > $altmpath/PKG_ADD
+    if test -z "`cat $altmpath/PKG_ADD`" ; then rm $altmpath/PKG_ADD; fi
+fi
+
+# grab the data files, skipping the CVS directory
+files=`echo $source/data/* | sed -e "s/[^ ]*CVS//"`
+if test -n "$files" -a "$files" != "$source/data/*" ; then
+    $INSTALL -d $mpath
+    $INSTALL_DATA $files $mpath
+fi
+
+# grab the executable files, skipping the CVS directory
+files=`echo $source/bin/* | sed -e "s/[^ ]*CVS//"`
+if test -n "$files" -a "$files" != "$source/bin/*" ; then
+    $INSTALL -d $xpath
+    $INSTALL_PROGRAM $files $xpath
+fi
+
+# grab the script files, skipping the CVS directory
+files=`echo $source/scripts/* | sed -e "s/[^ ]*CVS//"`
+if test -n "$files" -a "$files" != "$source/scripts/*" ; then
+    $INSTALL -d $xpath
+    $INSTALL_SCRIPT $files $xpath
+fi
+
--- a/main/linear-algebra/thfm.m	Sun Aug 20 14:44:45 2006 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,138 +0,0 @@
-%USAGE  y = thfm ( x, MODE )
-%
-%       trigonometric/hyperbolic functions of square matrix x
-%
-%MODE	cos   sin   tan   sec   csc   cot
-%	cosh  sinh  tanh  sech  csch  coth
-%       acos  asin  atan  asec  acsc  acot
-%       acosh asinh atanh asech acsch acoth
-%       sqrt  log   exp
-%
-%NOTE		--- IMPORTANT ---
-%	This algorithm does  NOT USE an eigensystem
-%	similarity transformation. It maps the MODE
-%	functions to functions of expm, logm and sqrtm, 
-%       which are known to be robust with respect to
-%	non-diagonalizable ('defective') x
-%
-%EXA	thfm( x ,'cos' )  calculates  matrix cosinus
-%	EVEN IF input matrix x IS NOT DIAGONALIZABLE
-%
-%ASSOC	expm, sqrtm, logm, funm
-% Copyright	(C) 2001 Rolf Fabian <fabian@tu-cottbus.de>
-% 010213
-%	published under current GNU GENERAL PUBLIC LICENSE
-
-% 2001-03-15 Paul Kienzle
-%     * extend with inverse functions and power functions
-%     * optimize handling of real input
-
-function y=thfm(x,M)
-				#% minimal arg check only
-  if	nargin~=2||~ischar(M)||ischar(x)	
-    usage ("y = thfm (x, MODE)");
-  endif
-
-  ## look for known functions of sqrt, log, exp
-  I = eye(size(x));
-  match = 1;
-  len =  length(M);
-  if	len==3
-    
-    if	M=='cos',  
-      if (isreal(x))     y = real( expm( i*x ) );
-      else               y = ( expm( i*x ) + expm( -i*x ) ) / 2;
-      endif
-      
-    elseif	M=='sin',
-      if (isreal(x))     y = imag( expm( i*x ) );
-      else               y = ( expm( i*x ) - expm( -i*x ) ) / (2*i);
-      endif
-      
-    elseif	M=='tan',
-      if (isreal(x))     t = expm( i*x );    y = imag(t)/real(t);
-      else     	         t = expm( -2*i*x ); y = -i*(I-t)/(I+t);
-      endif
-      
-    elseif	M=='cot',		% == cos/sin
-      if (isreal(x))     t = expm( i*x );    y = real(t)/imag(t);
-      else	         t = expm( -2*i*x ); y = i*(I+t)/(I-t);
-      endif
-      
-    elseif	M=='sec',  
-      if (isreal(x))     y = inv( real(expm(i*x)) );
-      else               y = inv( expm(i*x)+expm(-i*x) )*2 ;
-      endif
-      
-    elseif	M=='csc',  
-      if (isreal(x))     y = inv( imag(expm(i*x)) );
-      else               y = inv( expm(i*x)-expm(-i*x) )*2i;
-      endif
-
-    elseif    M=='log',  y = logm(x);
-      
-    elseif    M=='exp',  y = expm(x);
-      
-    else match = 0;
-
-    endif
-    
-  elseif	len==4
-    
-    if      M=='cosh',   y = ( expm(x)+expm(-x) )/2;
-      
-    elseif  M=='sinh',   y = ( expm(x)-expm(-x) )/2;
-      
-    elseif  M=='tanh'    t = expm( -2*x ); y = (I - t)/(I + t);
-     
-    elseif  M=='coth', 	 t = expm( -2*x ); y = (I + t)/(I - t);
-      
-    elseif  M=='sech',   y = 2 * inv( expm(x) + expm(-x) );
-      
-    elseif  M=='csch',   y = 2 * inv( expm(x) - expm(-x) );
-      
-    elseif  M=='asin',   y = -i * logm( i*x + sqrtm(I - x*x) );
-      
-    elseif  M=='acos',   y =  i * logm( x - i*sqrtm(I - x*x) );
-
-    elseif  M=='atan',   y = -i/2 * logm( (I + i*x)/(I - i*x) );
-
-    elseif  M=='acot',   y =  i/2 * logm( (I + i*x)/(i*x - I) );
-
-    elseif  M=='asec',   y = i * logm( ( I - sqrtm(I - x*x) ) / x );
-
-    elseif  M=='acsc',   y = -i * logm( i*( I + sqrtm(I - x*x) ) / x );
-
-    elseif  M=='sqrt',   y = sqrtm(x);
-
-    else match = 0;
-
-    end
-
-  elseif   len==5
-
-    if      M=='asinh',  y = logm( x + sqrtm (x*x + I) );
-      
-    elseif  M=='acosh',  y = logm( x + sqrtm (x*x - I) );
-      
-    elseif  M=='atanh',  y = logm( (I + x)/(I - x) ) / 2;
-      
-    elseif  M=='acoth',  y = logm( (I + x)/(x - I) ) / 2;
-
-    elseif  M=='asech',  y = logm( (I + sqrtm (I - x*x)) / x );
-
-    elseif  M=='acsch',  y = logm( (I + sqrtm (I + x*x)) / x );
-
-    else match = 0;
-
-    endif
-
-  else match = 0;
-
-  endif
-
-  ## if no known function found, use generic solver
-  if (match == 0)
-    y = funm( x, M );
-  endif
-endfunction