Mercurial > forge
changeset 2385:58386d13b8f1 octave-forge
Changed the directory structure of linear-algebra to match the package system
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