Mercurial > forge
changeset 6644:934b6c24e196 octave-forge
Remove the ARPACK package as it has been moved to core Octave
| author | hauberg |
|---|---|
| date | Fri, 05 Feb 2010 17:13:49 +0000 |
| parents | 9a49178a1115 |
| children | ea9af324d95d |
| files | nonfree/arpack/COPYING nonfree/arpack/DESCRIPTION nonfree/arpack/INDEX nonfree/arpack/doc/README nonfree/arpack/doc/Sorensen.eml nonfree/arpack/doc/patch nonfree/arpack/doc/patch.changelog nonfree/arpack/inst/svds.m nonfree/arpack/src/.svnignore nonfree/arpack/src/Makeconf.in nonfree/arpack/src/Makefile nonfree/arpack/src/autogen.sh nonfree/arpack/src/configure.base nonfree/arpack/src/eigs-base.cc nonfree/arpack/src/eigs.cc |
| diffstat | 15 files changed, 0 insertions(+), 6712 deletions(-) [+] |
line wrap: on
line diff
--- a/nonfree/arpack/COPYING Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,337 +0,0 @@ - GNU GENERAL PUBLIC LICENSE - Version 2, June 1991 - - Copyright (C) 1989, 1991 Free Software Foundation, Inc. <http://fsf.org/> - 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) <year> <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, see <http://www.gnu.org/licenses/>. - -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) year 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/nonfree/arpack/DESCRIPTION Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -Name: Arpack -Version: 1.0.8 -Date: 2009-05-03 -Author: David Bateman -Maintainer: David Bateman -Title: Arpack. -Description: Octave bindings to ARPACK, including the eigs and svds function. -Depends: octave (> 2.9.9) -Autoload: yes -License: GPL version 2 or later. Note ARPACK redistristibution - restriction discussed in the documentation -Url: http://octave.sf.net http://www.caam.rice.edu/software/ARPACK
--- a/nonfree/arpack/INDEX Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -sparse >> Sparse matrix support -Linear algebra - eigs - svds
--- a/nonfree/arpack/doc/README Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,54 +0,0 @@ -LICENSE: -------- - -Note that the ARPACK license includes the clause - -<quote> -Written notification is provided to the developers of intent to use this -software. Also, we ask that use of ARPACK is properly cited in any -resulting publications or software documentation. -</quote> - -This clause is difficult to interpret and so clarification has been -requested from the authors. Initial clarification can be found in the -e-mail Sorenson.eml. However, this is not considered to go far enough -to allow inclusion in Fedora. See the bug report - -https://bugzilla.redhat.com/bugzilla/show_bug.cgi?id=214967 - -for more information. - -However, the e-mail does make it clear that there is no restrictions -on the distribution of this binding of ARPACK to Octave, and no -restrictions on its use. Furthermore, I take the clarification to mean -that modifications of these Octave bindings are allowed without -further notification to the authors. However, if ARPACK is to be used -with another package, even if it uses the octave code as a basis, then -further notification must be obtained. - -DOWNLOAD: --------- - -The main webpage for ARPACK is http://www.caam.rice.edu/software/ARPACK/ -where two files can be obtained, these being - -http://www.caam.rice.edu/software/ARPACK/SRC/arpack96.tar.gz -http://www.caam.rice.edu/software/ARPACK/SRC/patch.tar.gz - -Note that the second file is not a patch, but rather another tar -archive that overwrites the existing files with replacements - -Chao Yang maintains an updated version of ARPACK, and is one of the -authors of ARPACK. However, I'm unsure if this version has any -official status. It can be found at - -http://www.cs.ucsb.edu/~viral/parpack96/patch.tar.gz - -PATCHES: -------- -Included here are a patch and a changelog file that were created -against the CVS of octave as of 22 Feb 2007. These might be used to -port this code to Octave once the above licensing question is -addressed. Note that the code eigs.cc, eigs-base.cc and svds.m must -also be copied to an appropriate place in the Octave source tree when -applying this patch
--- a/nonfree/arpack/doc/Sorensen.eml Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,102 +0,0 @@ -X-Mozilla-Status: 0011 -X-Mozilla-Status2: 00000000 -Received: from zuk35exm65.ds.mot.com ([10.178.1.44]) by zuk35exm62.ds.mot.com with Microsoft SMTPSVC(6.0.3790.2709); - Tue, 13 Feb 2007 19:53:15 +0000 -Received: from az33exr02.mot.com ([10.64.251.232]) by zuk35exm65.ds.mot.com with Microsoft SMTPSVC(6.0.3790.2709); - Tue, 13 Feb 2007 19:53:14 +0000 -Received: from motgate3.mot.com (motgate3.mot.com [144.189.100.103]) - by az33exr02.mot.com (8.13.1/8.13.0) with ESMTP id l1DJrDDS020683 - for <David.Bateman@motorola.com>; Tue, 13 Feb 2007 13:53:13 -0600 (CST) -Received: from mail128.messagelabs.com (mail128.messagelabs.com [216.82.250.131]) - by motgate3.mot.com (8.12.11/Motorola) with SMTP id l1DJr8gl023906 - for <David.Bateman@motorola.com>; Tue, 13 Feb 2007 12:53:11 -0700 (MST) -X-VirusChecked: Checked -X-Env-Sender: sorensen@rice.edu -X-Msg-Ref: server-6.tower-128.messagelabs.com!1171396385!13383447!1 -X-StarScan-Version: 5.5.10.7.1; banners=-,-,- -X-Originating-IP: [128.42.17.10] -X-SpamReason: No, hits=0.0 required=7.0 tests= -Received: (qmail 3872 invoked from network); 13 Feb 2007 19:53:05 -0000 -Received: from caam.rice.edu (HELO caam.rice.edu) (128.42.17.10) - by server-6.tower-128.messagelabs.com with SMTP; 13 Feb 2007 19:53:05 -0000 -Received: from localhost (localhost [127.0.0.1]) - by caam.rice.edu (Postfix) with ESMTP id 64341153A7 - for <David.Bateman@motorola.com>; Tue, 13 Feb 2007 13:53:04 -0600 (CST) -Received: from caam.rice.edu ([127.0.0.1]) - by localhost (caam.rice.edu [127.0.0.1]) (amavisd-new, port 10024) with LMTP - id 23777-01-16 for <David.Bateman@motorola.com>; - Tue, 13 Feb 2007 13:52:59 -0600 (CST) -Received: from [128.42.21.177] (sorensenl400.caam.rice.edu [128.42.21.177]) - by caam.rice.edu (Postfix) with ESMTP id 494E81539F - for <David.Bateman@motorola.com>; Tue, 13 Feb 2007 13:52:59 -0600 (CST) -Message-ID: <45D2171B.8030109@rice.edu> -Date: Tue, 13 Feb 2007 13:52:59 -0600 -From: Dan Sorensen<sorensen@rice.edu> -User-Agent: Thunderbird 1.5.0.9 (Windows/20061207) -MIME-Version: 1.0 -To: David Bateman<David.Bateman@motorola.com> -Subject: Re: ARPACK License Question -References: <457EE5B3.70402@ieee.org> <20070105114426.GI4860@neu.nirvana> <45B8CB2F.9030904@motorola.com> -In-Reply-To: <45B8CB2F.9030904@motorola.com> -Content-Type: text/plain; charset=ISO-8859-1; format=flowed -Content-Transfer-Encoding: 7bit -X-Virus-Scanned: by amavis-2.2.1 at caam.rice.edu -Return-Path: sorensen@rice.edu -X-OriginalArrivalTime: 13 Feb 2007 19:53:14.0831 (UTC) FILETIME=[994B7DF0:01C74FA8] - -Dear Mr. Bateman - -I apologize for not responding to this previously. - -The clarification we discussed is the following - - -The clause in the license statement that states - ->>Written notification is provided to the developers of intent to use this ->> software. Also, we ask that use of ARPACK is properly cited in any ->> resulting publications or software documentation. - -has the following intension in your case. - -We are asking for acknowledgment in FEDORA that ARPACK is -the software that underlies what corresponds to the "eigs" command. -There is no intention to pass on a requirement of notification of use -from users of FEDORA. - -This is the understanding we have with MATLAB for example. - -If the above note or a slight modification of it is not acceptable -for the purposes of using ARPACK in FEDORA, I will have to refer -you to the tech transfer department of Rice University as I explained -during our phone conversation. - -Once again my apologies for the delay and I thank you for your -interest in ARPACK. - -Best Regards -Dan Sorensen - - - - - - - -David Bateman wrote: -> Dear Professor Sorensen, -> -> Perhaps you have not yet seen the e-mail below, and so I draw it to your -> attention. Can you please examine the request to modify the license of -> ARPACK in this mail belong to allow its inclusion in FEDORA and other -> similar open source linux distributions? -> -> As the author of the eigs function for Octave (www.octave.org) that uses -> ARPACK for its functionality, I'd hate to see my work not included in -> Octave due to this question not being resolved. -> -> Best Regards -> David -> -> -
--- a/nonfree/arpack/doc/patch Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,264 +0,0 @@ -*** ./Makeconf.in.orig41 2007-02-21 13:19:49.648245987 +0100 ---- ./Makeconf.in 2007-02-21 13:39:38.661764244 +0100 -*************** -*** 221,226 **** ---- 221,227 ---- - CCOLAMD_LIBS = @CCOLAMD_LIBS@ - CHOLMOD_LIBS = @CHOLMOD_LIBS@ - CXSPARSE_LIBS = @CXSPARSE_LIBS@ -+ ARPACK_LIBS = @ARPACK_LIBS@ - LIBS = @LIBS@ - - USE_64_BIT_IDX_T = @USE_64_BIT_IDX_T@ -*** ./doc/interpreter/sparse.txi.orig41 2007-02-21 13:33:25.269502039 +0100 ---- ./doc/interpreter/sparse.txi 2007-02-21 13:38:33.886013421 +0100 -*************** -*** 421,430 **** - @dfn{csymamd}, @dfn{dmperm}, @dfn{symamd}, @dfn{randperm}, (symrcm) - - @item Linear algebra: -! @dfn{matrix\_type}, @dfn{spchol}, @dfn{cpcholinv}, -! @dfn{spchol2inv}, @dfn{spdet}, @dfn{spinv}, @dfn{spkron}, -! @dfn{splchol}, @dfn{splu}, @dfn{spqr}, (condest, eigs, normest, -! sprank, svds, spaugment) - - @item Iterative techniques: - @dfn{luinc}, @dfn{pcg}, @dfn{pcr}, (bicg, bicgstab, cholinc, cgs, ---- 421,430 ---- - @dfn{csymamd}, @dfn{dmperm}, @dfn{symamd}, @dfn{randperm}, (symrcm) - - @item Linear algebra: -! @dfn{eigs}, @dfn{matrix\_type}, @dfn{normest}, @dfn{spchol}, -! @dfn{spcholinv}, @dfn{spchol2inv}, @dfn{spdet}, @dfn{spinv}, @dfn{spkron}, -! @dfn{splchol}, @dfn{splu}, @dfn{spqr}, @dfn{sprank}, @dfn{svds}, -! (condest, spaugment) - - @item Iterative techniques: - @dfn{luinc}, @dfn{pcg}, @dfn{pcr}, (bicg, bicgstab, cholinc, cgs, -*************** -*** 1434,1440 **** - @item condest - @emph{Not implemented} - @item eigs -! @emph{Not implemented} - @item @ref{matrix_type} - Identify the matrix type or mark a matrix as a particular type. - @item @ref{normest} ---- 1434,1441 ---- - @item condest - @emph{Not implemented} - @item eigs -! Calculate a limited number of eigenvalues and eigenvectors of @var{a}, -! based on a selection criteria. - @item @ref{matrix_type} - Identify the matrix type or mark a matrix as a particular type. - @item @ref{normest} -*************** -*** 1462,1468 **** - @item @ref{sprank} - Calculates the structural rank of a sparse matrix @var{s}. - @item svds -! @emph{Not implemented} - @end table - @subsubsection Iterative techniques - @table @asis ---- 1463,1469 ---- - @item @ref{sprank} - Calculates the structural rank of a sparse matrix @var{s}. - @item svds -! Find a few singular values of the matrix @var{a}. - @end table - @subsubsection Iterative techniques - @table @asis -*************** -*** 1537,1542 **** ---- 1538,1545 ---- - sparser Cholesky factor than S. - * dmperm:: Perfrom a Deulmage-Mendelsohn permutation on the sparse - matrix S. -+ * eigs:: Calculate a limited number of eigenvalues and eigenvectors -+ of @var{a}, based on a selection criteria. - * etree:: Returns the elimination tree for the matrix S. - * etreeplot:: Plots the elimination tree of the matrix @var{s} or - @code{@var{s}+@var{s}'} if @var{s} in non-symmetric. -*************** -*** 1608,1613 **** ---- 1611,1617 ---- - * spsum:: Sum of elements along dimension DIM. - * spsumsq:: Sum of squares of elements along dimension DIM. - * spy:: Plot the sparsity pattern of the sparse matrix X -+ * svds:: Find a few singular values of the matrix @var{a}. - * symamd:: For a symmetric positive definite matrix S, returns the - permutation vector p such that `S (P, P)' tends to have a - sparser Cholesky factor than S. -*************** -*** 1636,1647 **** - - @DOCSTRING(csymamd) - -! @node dmperm, etree, csymamd, Function Reference - @subsubsection dmperm - - @DOCSTRING(dmperm) - -! @node etree, etreeplot, dmperm, Function Reference - @subsubsection etree - - @DOCSTRING(etree) ---- 1640,1656 ---- - - @DOCSTRING(csymamd) - -! @node dmperm, eigs, csymamd, Function Reference - @subsubsection dmperm - - @DOCSTRING(dmperm) - -! @node eigs, etree, dmperm, Function Reference -! @subsubsection eigs -! -! @DOCSTRING(eigs) -! -! @node etree, etreeplot, eigs, Function Reference - @subsubsection etree - - @DOCSTRING(etree) -*************** -*** 1867,1878 **** - - @DOCSTRING(spsumsq) - -! @node spy, symamd, spsumsq, Function Reference - @subsubsection spy - - @DOCSTRING(spy) - -! @node symamd, symbfact, spy, Function Reference - @subsubsection symamd - - @DOCSTRING(symamd) ---- 1876,1892 ---- - - @DOCSTRING(spsumsq) - -! @node spy, svds, spsumsq, Function Reference - @subsubsection spy - - @DOCSTRING(spy) - -! @node svds, symamd, spy, Function Reference -! @subsubsection svds -! -! @DOCSTRING(svds) -! -! @node symamd, symbfact, svds, Function Reference - @subsubsection symamd - - @DOCSTRING(symamd) -*** ./src/Makefile.in.orig41 2007-02-21 13:19:49.646246085 +0100 ---- ./src/Makefile.in 2007-02-21 13:39:38.652764696 +0100 -*************** -*** 49,56 **** - - DLD_XSRC := balance.cc besselj.cc betainc.cc cellfun.cc chol.cc \ - ccolamd.cc colamd.cc colloc.cc conv2.cc daspk.cc dasrt.cc \ -! dassl.cc det.cc dispatch.cc eig.cc expm.cc fft.cc fft2.cc \ -! fftn.cc fftw.cc filter.cc find.cc fsolve.cc \ - gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \ - givens.cc hess.cc interpn.cc inv.cc kron.cc lpsolve.cc lsode.cc \ - lu.cc luinc.cc matrix_type.cc minmax.cc pinv.cc qr.cc \ ---- 49,56 ---- - - DLD_XSRC := balance.cc besselj.cc betainc.cc cellfun.cc chol.cc \ - ccolamd.cc colamd.cc colloc.cc conv2.cc daspk.cc dasrt.cc \ -! dassl.cc det.cc dispatch.cc eig.cc eigs.cc expm.cc fft.cc \ -! fft2.cc fftn.cc fftw.cc filter.cc find.cc fsolve.cc \ - gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \ - givens.cc hess.cc interpn.cc inv.cc kron.cc lpsolve.cc lsode.cc \ - lu.cc luinc.cc matrix_type.cc minmax.cc pinv.cc qr.cc \ -*************** -*** 577,582 **** ---- 577,584 ---- - - ifeq ($(ENABLE_DYNAMIC_LINKING), true) - ifdef CXXPICFLAG -+ eigs.oct : pic/eigs.o octave$(EXEEXT) -+ $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(ARPACK_LIBS) - regexp.oct : pic/regexp.o octave$(EXEEXT) - $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(REGEX_LIBS) - urlwrite.oct : pic/urlwrite.o octave$(EXEEXT) -*************** -*** 584,589 **** ---- 586,593 ---- - __glpk__.oct : pic/__glpk__.o octave$(EXEEXT) - $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(GLPK_LIBS) - else -+ eigs.oct : eigs.o octave$(EXEEXT) -+ $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(ARPACK_LIBS) - regexp.oct : regexp.o octave$(EXEEXT) - $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(REGEX_LIBS) - urlwrite.oct : urlwrite.o octave$(EXEEXT) -*** ./configure.in.orig41 2007-02-21 13:19:49.648245987 +0100 ---- ./configure.in 2007-02-21 13:39:38.659764345 +0100 -*************** -*** 1001,1006 **** ---- 1001,1035 ---- - AC_MSG_WARN($warn_cxsparse) - fi - -+ WITH_ARPACK=true -+ AC_ARG_WITH(arpack, -+ [AS_HELP_STRING([--without-arpack], [don't use arpack])], -+ with_arpack=$withval, with_arpack=yes) -+ -+ arpack_lib= -+ if test "$with_arpack" = yes; then -+ arpack_lib="arpack" -+ elif test "$with_arpack" != no; then -+ arpack_lib="$with_arpack" -+ fi -+ -+ ARPACK_LIBS= -+ AC_SUBST(ARPACK_LIBS) -+ WITH_ARPACK=false -+ if test -n "$arpack_lib"; then -+ AC_CHECK_LIB($arpack_lib, F77_FUNC(dseupd,DSEUPD), [ -+ WITH_ARPACK=true -+ ARPACK_LIBS="-l$arpack_lib" -+ FLIBS="$ARPACK_LIBS $FLIBS" -+ AC_DEFINE(HAVE_ARPACK, 1, [Define if ARPACK is available.])], , -+ $BLAS_LIBS $FLIBS) -+ fi -+ if test $WITH_ARPACK = no; then -+ warn_arpack="arpack not found. This will result in a lack of the eigs function." -+ AC_MSG_WARN($warn_arpack) -+ fi -+ -+ - ### Handle shared library options. - - ### Enable creation of static libraries. -*************** -*** 2019,2024 **** ---- 2048,2054 ---- - CCOLAMD libraries: $CCOLAMD_LIBS - CHOLMOD libraries: $CHOLMOD_LIBS - CXSPARSE libraries: $CXSPARSE_LIBS -+ ARPACK libraries: $ARPACK_LIBS - HDF5 libraries: $HDF5_LIBS - CURL libraries: $CURL_LIBS - REGEX libraries: $REGEX_LIBS -*************** -*** 2126,2131 **** ---- 2156,2166 ---- - warn_msg_printed=true - fi - -+ if test -n "$warn_arpack"; then -+ AC_MSG_WARN($warn_arpack) -+ warn_msg_printed=true -+ fi -+ - if test -n "$warn_fftw"; then - AC_MSG_WARN($warn_fftw) - warn_msg_printed=true
--- a/nonfree/arpack/doc/patch.changelog Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,24 +0,0 @@ -2007-02-21 David Bateman <dbateman@free.fr> - - * interpreter/sparse.txi: Document the svds and eigs functions - -2007-02-21 David Bateman <dbateman@free.fr> - - * configure.in (ARPACK_LIBS): Add test for ARPACK. Warn if it - isn't found. Set HAVE_ARPACK and ARPACK_LIBS. - * Makeconf.in (ARPACK_LIBS): New variable. - -2007-02-21 David Bateman <dbateman@free.fr> - - * Makefile.in: Add eigs-base.cc to TEMPLATE_SRC. - * eigs-base.cc: New template class file for ARPACk functions. - -2007-02-21 David Bateman <dbateman@free.fr> - - * Makefile.in: Add specific build for eigs.oct to include -larpack. - (DLD_XSRC): Add eigs.cc. - * DLD-FUNCTIONS/eigs.cc: New file. - -2007-02-21 David Bateman <dbateman@free.fr> - - * sparse/svds.m: New file.
--- a/nonfree/arpack/inst/svds.m Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,231 +0,0 @@ -## Copyright (C) 2006 David Bateman -## -## 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, see <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn {Function File} {@var{s} =} svds (@var{a}) -## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k}) -## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k}, @var{sigma}) -## @deftypefnx {Function File} {@var{s} =} svds (@var{a}, @var{k}, @var{sigma}, @var{opts}) -## @deftypefnx {Function File} {[@var{u}, @var{s}, @var{v}, @var{flag}] =} svds (@dots{}) -## -## Find a few singular values of the matrix @var{a}. The singular values -## are calculated using -## -## @example -## @group -## [@var{m}, @var{n}] = size(@var{a}) -## @var{s} = eigs([sparse(@var{m}, @var{m}), @var{a}; @dots{} -## @var{a}', sparse(@var{n}, @var{n})]) -## @end group -## @end example -## -## The eigenvalues returned by @code{eigs} correspond to the singular -## values of @var{a}. The number of singular values to calculate is given -## by @var{k}, whose default value is 6. -## -## The argument @var{sigma} can be used to specify which singular values -## to find. @var{sigma} can be either the string 'L', the default, in -## which case the largest singular values of @var{a} are found. Otherwise -## @var{sigma} should be a real scalar, in which case the singular values -## closest to @var{sigma} are found. Note that for relatively small values -## of @var{sigma}, there is the chance that the requested number of singular -## values are not returned. In that case @var{sigma} should be increased. -## -## If @var{opts} is given, then it is a structure that defines options -## that @code{svds} will pass to @var{eigs}. The possible fields of this -## structure are therefore determined by @code{eigs}. By default three -## fields of this structure are set by @code{svds}. -## -## @table @code -## @item tol -## The required convergence tolerance for the singular values. @code{eigs} -## is passed @var{tol} divided by @code{sqrt(2)}. The default value is -## 1e-10. -## -## @item maxit -## The maximum number of iterations. The defaut is 300. -## -## @item disp -## The level of diagnostic printout. If @code{disp} is 0 then there is no -## printout. The default value is 0. -## @end table -## -## If more than one output argument is given, then @code{svds} also -## calculates the left and right singular vectors of @var{a}. @var{flag} -## is used to signal the convergence of @code{svds}. If @code{svds} -## converges to the desired tolerance, then @var{flag} given by -## -## @example -## @group -## norm (@var{a} * @var{v} - @var{u} * @var{s}, 1) <= @dots{} -## @var{tol} * norm (@var{a}, 1) -## @end group -## @end example -## -## will be zero. -## @end deftypefn -## @seealso{eigs} - -function [u, s, v, flag] = svds (a, k, sigma, opts) - - if (nargin < 1 || nargin > 4) - error ("Incorrect number of arguments"); - endif - - if (nargin < 4) - opts.tol = 1e-10 / sqrt(2); - opts.disp = 0; - opts.maxit = 300; - else - if (!isstruct(opts)) - error("opts must be a structure"); - endif - if (!isfield(opts,"tol")) - opts.tol = 1e-10 / sqrt(2); - endif - endif - - if (nargin < 3 || strcmp(sigma,"L")) - if (isreal(a)) - sigma = "LA"; - else - sigma = "LR"; - endif - elseif (isscalar(sigma) && isreal(sigma)) - if ((sigma < 0)) - error ("sigma must be a positive real value"); - endif - else - error ("sigma must be a positive real value or the string 'L'"); - endif - - maxA = max(max(abs(a))); - if (maxA == 0) - u = eye(m, k); - s = zeros(k, k); - v = eye(n, k); - else - [m, n] = size(a); - if (nargin < 2) - k = min([6, m, n]); - else - k = min([k, m, n]); - endif - - ## Scale everything by the 1-norm to make things more stable. - B = a / maxA; - Bopts = opts; - Bopts.tol = opts.tol / maxA; - Bsigma = sigma; - if (!ischar(Bsigma)) - Bsigma = Bsigma / maxA; - endif - - if (!ischar(Bsigma) && Bsigma == 0) - ## The eigenvalues returns by eigs are symmetric about 0. As we - ## are only interested in the positive eigenvalues, we have to - ## double k. If sigma is smaller than the smallest singular value - ## this can also be an issue. However, we'd like to avoid double - ## k for all scalar value of sigma... - [V, s, flag] = eigs ([sparse(m,m), B; B', sparse(n,n)], - 2 * k, Bsigma, Bopts); - else - [V, s, flag] = eigs ([sparse(m,m), B; B', sparse(n,n)], - k, Bsigma, Bopts); - endif - s = diag(s); - - if (ischar(sigma)) - norma = max(s); - else - norma = normest(a); - endif - V = sqrt(2) * V; - u = V(1:m,:); - v = V(m+1:end,:); - - ## We wish to exclude all eigenvalues that are less than zero as these - ## are artifacts of the way the matrix passed to eigs is formed. There - ## is also the possibility that the value of sigma chosen is exactly - ## a singular value, and in that case we're dead!! So have to rely on - ## the warning from eigs. We exclude the singular values which are - ## less than or equal to zero to within some tolerance scaled by the - ## norm since if we don't we might end up with too many singular - ## values. What is appropriate for the tolerance? - tol = norma * opts.tol; - ind = find(s > tol); - if (length(ind) < k) - ## Find the zero eigenvalues of B, Ignore the eigenvalues that are - ## nominally negative. - zind = find(abs(s) <= tol); - p = min(length(zind), k-length(ind)); - ind = [ind;zind(1:p)]; - elseif (length(ind) > k) - ind = ind(1:k); - endif - u = u(:,ind); - s = s(ind); - v = v(:,ind); - - if (length(s) < k) - warning("returning fewer singular values than requested."); - if (!ischar(sigma)) - warning("try increasing the value of sigma"); - endif - endif - - s = s * maxA; - endif - - if (nargout < 2) - u = s; - else - s = diag(s); - if (nargout > 3) - flag = norm(a*v - u*s, 1) > sqrt(2) * opts.tol * norm(a, 1); - endif - endif -endfunction - -%!shared n, k, a, u, s, v, opts -%! n = 100; -%! k = 7; -%! a = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),0.4*n*ones(1,n),ones(1,n-2)]); -%! %%a = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]); -%! [u,s,v] = svd(full(a)); -%! s = diag(s); -%! [dum, idx] = sort(abs(s)); -%! s = s(idx); -%! u = u(:,idx); -%! v = v(:,idx); -%! randn('state',42) -%!test -%! [u2,s2,v2,flag] = svds(a,k); -%! s2 = diag(s2); -%! assert(flag,!1); -%! assert(s(end:-1:end-k+1), s2, 1e-10); -%!test -%! [u2,s2,v2,flag] = svds(a,k,0); -%! s2 = diag(s2); -%! assert(flag,!1); -%! assert(s(k:-1:1), s2, 1e-10); -%!test -%! idx = floor(n/2); -%! % Don't put sigma right on a singular value or there are convergence -%! sigma = 0.99*s(idx) + 0.01*s(idx+1); -%! [u2,s2,v2,flag] = svds(a,k,sigma); -%! s2 = diag(s2); -%! assert(flag,!1); -%! assert(s((idx+floor(k/2)):-1:(idx-floor(k/2))), s2, 1e-10);
--- a/nonfree/arpack/src/.svnignore Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,2 +0,0 @@ -configure -autom4te.cache
--- a/nonfree/arpack/src/Makeconf.in Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,70 +0,0 @@ - -## 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@ -datarootdir = @datarootdir@ -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@ - -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@ - -ver = @ver@ -MPATH = @mpath@ -OPATH = @opath@ -XPATH = @xpath@ -ALTMPATH = @altmpath@ -ALTOPATH = @altopath@ - -HAVE_ARPACK=@HAVE_ARPACK@ -ARPACK_LIBS=@ARPACK_LIBS@ - -%.o: %.c ; $(MKOCTFILE) -c $< -%.o: %.f ; $(MKOCTFILE) -c $< -%.o: %.cc ; $(MKOCTFILE) -c $< -%.oct: %.cc ; $(MKOCTFILE) $<
--- a/nonfree/arpack/src/Makefile Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,12 +0,0 @@ -sinclude Makeconf - -ifeq ($(HAVE_ARPACK),yes) - PROGS:=eigs.oct -endif - -all: $(PROGS) - -eigs.oct: eigs.cc - $(MKOCTFILE) -DHAVE_CONFIG_H -DHAVE_ARPACK $< $(ARPACK_LIBS) - -clean: ; -$(RM) *.o core octave-core *.oct *~
--- a/nonfree/arpack/src/autogen.sh Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,27 +0,0 @@ -#! /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 -cat <<EOF >> configure.in - AC_OUTPUT(\$CONFIGURE_OUTPUTS) - dnl XXX FIXME XXX chmod is not in autoconf's list of portable functions - - 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 configure.in > configure.tmp -if [ diff configure.tmp configure > /dev/null 2>&1 ]; then - rm -f configure.tmp; -else - mv -f configure.tmp configure - chmod 0755 configure -fi -rm -f configure.in
--- a/nonfree/arpack/src/configure.base Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,350 +0,0 @@ -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,extra_libs) -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 $5" -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 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_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 - -AC_ARG_WITH(arpack, - [AS_HELP_STRING([--without-arpack], [don't use arpack])], - with_arpack=$withval, with_arpack=yes) - -arpack_lib= -if test "$with_arpack" = yes; then - arpack_lib="arpack" -elif test "$with_arpack" != no; then - arpack_lib="$with_arpack" -fi - -ARPACK_LIBS= -HAVE_ARPACK=no -AC_SUBST(ARPACK_LIBS) -AC_SUBST(HAVE_ARPACK) -WITH_ARPACK=no -if test -n "$arpack_lib"; then - TRY_MKOCTFILE([for libarpack],[char F77_FUNC(dseupd,DSEUPD) ();],[ - WITH_ARPACK=yes - ARPACK_LIBS="-l$arpack_lib" - HAVE_ARPACK=yes], ,-l$arpack_lib) -fi - -if test $WITH_ARPACK = no; then - AC_MSG_WARN("arpack not found. This essentially makes this package useless.") - ARPACKSTATUS="libarpack not found" -else - ARPACKSTATUS="Using ARPACK libraries $ARPACK_LIBS" -fi - -CONFIGURE_OUTPUTS="Makeconf" -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 - arpack: $ARPACKSTATUS"
--- a/nonfree/arpack/src/eigs-base.cc Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,3760 +0,0 @@ -/* - -Copyright (C) 2005 David Bateman - -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 software 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 software; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <cfloat> -#include <cmath> -#include <vector> -#include <iostream> - -#include "f77-fcn.h" -#include "quit.h" -#include "SparsedbleLU.h" -#include "SparseCmplxLU.h" -#include "dSparse.h" -#include "CSparse.h" -#include "MatrixType.h" -#include "SparsedbleCHOL.h" -#include "SparseCmplxCHOL.h" -#include "oct-rand.h" -#include "dbleCHOL.h" -#include "CmplxCHOL.h" -#include "dbleLU.h" -#include "CmplxLU.h" - -#ifdef HAVE_ARPACK -typedef ColumnVector (*EigsFunc) (const ColumnVector &x, int &eigs_error); -typedef ComplexColumnVector (*EigsComplexFunc) - (const ComplexColumnVector &x, int &eigs_error); - -// Arpack and blas fortran functions we call. -extern "C" -{ - F77_RET_T - F77_FUNC (dsaupd, DSAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const double&, - double*, const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type*, - octave_idx_type*, double*, double*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (dseupd, DSEUPD) (const int&, F77_CONST_CHAR_ARG_DECL, - int*, double*, double*, - const octave_idx_type&, const double&, - F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, - const double&, double*, const octave_idx_type&, - double*, const octave_idx_type&, octave_idx_type*, - octave_idx_type*, double*, double*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (dnaupd, DNAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - octave_idx_type&, const double&, - double*, const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type*, - octave_idx_type*, double*, double*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (dneupd, DNEUPD) (const int&, F77_CONST_CHAR_ARG_DECL, - int*, double*, double*, - double*, const octave_idx_type&, const double&, - const double&, double*, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - octave_idx_type&, const double&, double*, - const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type*, - octave_idx_type*, double*, double*, - const octave_idx_type&, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (znaupd, ZNAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const double&, - Complex*, const octave_idx_type&, Complex*, - const octave_idx_type&, octave_idx_type*, - octave_idx_type*, Complex*, Complex*, - const octave_idx_type&, double *, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (zneupd, ZNEUPD) (const int&, F77_CONST_CHAR_ARG_DECL, - int*, Complex*, Complex*, - const octave_idx_type&, const Complex&, - Complex*, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const double&, - Complex*, const octave_idx_type&, Complex*, - const octave_idx_type&, octave_idx_type*, - octave_idx_type*, Complex*, Complex*, - const octave_idx_type&, double *, octave_idx_type& - F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL - F77_CHAR_ARG_LEN_DECL); - - F77_RET_T - F77_FUNC (dgemv, DGEMV) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, const double&, - const double*, const octave_idx_type&, const double*, - const octave_idx_type&, const double&, double*, - const octave_idx_type& - F77_CHAR_ARG_LEN_DECL); - - - F77_RET_T - F77_FUNC (zgemv, ZGEMV) (F77_CONST_CHAR_ARG_DECL, - const octave_idx_type&, const octave_idx_type&, const Complex&, - const Complex*, const octave_idx_type&, const Complex*, - const octave_idx_type&, const Complex&, Complex*, const octave_idx_type& - F77_CHAR_ARG_LEN_DECL); - -} - - -#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) -static octave_idx_type -lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&); - -static octave_idx_type -lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&, - ComplexMatrix&); - -static octave_idx_type -lusolve (const Matrix&, const Matrix&, Matrix&); - -static octave_idx_type -lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&); - -static ComplexMatrix -ltsolve (const SparseComplexMatrix&, const ColumnVector&, - const ComplexMatrix&); - -static Matrix -ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&,); - -static ComplexMatrix -ltsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&); - -static Matrix -ltsolve (const Matrix&, const ColumnVector&, const Matrix&,); - -static ComplexMatrix -utsolve (const SparseComplexMatrix&, const ColumnVector&, const ComplexMatrix&); - -static Matrix -utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&); - -static ComplexMatrix -utsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&); - -static Matrix -utsolve (const Matrix&, const ColumnVector&, const Matrix&); - -#endif - -template <class M, class SM> -static octave_idx_type -lusolve (const SM& L, const SM& U, M& m) -{ - octave_idx_type err = 0; - double rcond; - MatrixType utyp (MatrixType::Upper); - - // Sparse L is lower triangular, Dense L is permuted lower triangular!!! - m = L.solve (m, err, rcond, 0); - if (err) - return err; - - m = U.solve (utyp, m, err, rcond, 0); - - return err; -} - -template <class SM, class M> -static M -ltsolve (const SM& L, const ColumnVector& Q, const M& m) -{ - octave_idx_type n = L.cols(); - octave_idx_type b_nc = m.cols(); - octave_idx_type err = 0; - double rcond; - MatrixType ltyp (MatrixType::Lower); - M tmp = L.solve (ltyp, m, err, rcond, 0); - M retval; - const double* qv = Q.fortran_vec(); - - if (!err) - { - retval.resize (n, b_nc); - for (octave_idx_type j = 0; j < b_nc; j++) - { - for (octave_idx_type i = 0; i < n; i++) - retval.elem(static_cast<octave_idx_type>(qv[i]), j) = - tmp.elem(i,j); - } - } - - return retval; -} - -template <class SM, class M> -static M -utsolve (const SM& U, const ColumnVector& Q, const M& m) -{ - octave_idx_type n = U.cols(); - octave_idx_type b_nc = m.cols(); - octave_idx_type err = 0; - double rcond; - MatrixType utyp (MatrixType::Upper); - - M retval (n, b_nc); - const double* qv = Q.fortran_vec(); - for (octave_idx_type j = 0; j < b_nc; j++) - { - for (octave_idx_type i = 0; i < n; i++) - retval.elem(i,j) = m.elem(static_cast<octave_idx_type>(qv[i]), j); - } - return U.solve (utyp, retval, err, rcond, 0); -} - -static bool -vector_product (const SparseMatrix& m, const double* x, double* y) -{ - octave_idx_type nc = m.cols (); - - for (octave_idx_type j = 0; j < nc; j++) - y[j] = 0.; - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) - y[m.ridx(i)] += m.data(i) * x[j]; - - return true; -} - -static bool -vector_product (const Matrix& m, const double *x, double *y) -{ - octave_idx_type nr = m.rows (); - octave_idx_type nc = m.cols (); - - F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1), - nr, nc, 1.0, m.data (), nr, - x, 1, 0.0, y, 1 - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable error in dgemv"); - return false; - } - else - return true; -} - -static bool -vector_product (const SparseComplexMatrix& m, const Complex* x, - Complex* y) -{ - octave_idx_type nc = m.cols (); - - for (octave_idx_type j = 0; j < nc; j++) - y[j] = 0.; - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) - y[m.ridx(i)] += m.data(i) * x[j]; - - return true; -} - -static bool -vector_product (const ComplexMatrix& m, const Complex *x, Complex *y) -{ - octave_idx_type nr = m.rows (); - octave_idx_type nc = m.cols (); - - F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 ("N", 1), - nr, nc, 1.0, m.data (), nr, - x, 1, 0.0, y, 1 - F77_CHAR_ARG_LEN (1))); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable error in zgemv"); - return false; - } - else - return true; -} - -static bool -make_cholb (Matrix& b, Matrix& bt, ColumnVector& permB) -{ - octave_idx_type info; - CHOL fact (b, info); - octave_idx_type n = b.cols(); - - if (info != 0) - return false; - else - { - bt = fact.chol_matrix (); - b = bt.transpose(); - permB = ColumnVector(n); - for (octave_idx_type i = 0; i < n; i++) - permB(i) = i; - return true; - } -} - -static bool -make_cholb (SparseMatrix& b, SparseMatrix& bt, ColumnVector& permB) -{ - octave_idx_type info; - SparseCHOL fact (b, info, false); - - if (fact.P() != 0) - return false; - else - { - b = fact.L(); - bt = b.transpose(); - permB = fact.perm() - 1.0; - return true; - } -} - -static bool -make_cholb (ComplexMatrix& b, ComplexMatrix& bt, ColumnVector& permB) -{ - octave_idx_type info; - ComplexCHOL fact (b, info); - octave_idx_type n = b.cols(); - - if (info != 0) - return false; - else - { - bt = fact.chol_matrix (); - b = bt.hermitian(); - permB = ColumnVector(n); - for (octave_idx_type i = 0; i < n; i++) - permB(i) = i; - return true; - } -} - -static bool -make_cholb (SparseComplexMatrix& b, SparseComplexMatrix& bt, - ColumnVector& permB) -{ - octave_idx_type info; - SparseComplexCHOL fact (b, info, false); - - if (fact.P() != 0) - return false; - else - { - b = fact.L(); - bt = b.hermitian(); - permB = fact.perm() - 1.0; - return true; - } -} - -static bool -LuAminusSigmaB (const SparseMatrix &m, const SparseMatrix &b, - bool cholB, const ColumnVector& permB, double sigma, - SparseMatrix &L, SparseMatrix &U, octave_idx_type *P, - octave_idx_type *Q) -{ - bool have_b = ! b.is_empty (); - octave_idx_type n = m.rows(); - - // Caclulate LU decomposition of 'A - sigma * B' - SparseMatrix AminusSigmaB (m); - - if (have_b) - { - if (cholB) - { - if (permB.length()) - { - SparseMatrix tmp(n,n,n); - for (octave_idx_type i = 0; i < n; i++) - { - tmp.xcidx(i) = i; - tmp.xridx(i) = - static_cast<octave_idx_type>(permB(i)); - tmp.xdata(i) = 1; - } - tmp.xcidx(n) = n; - - AminusSigmaB = AminusSigmaB - sigma * tmp * - b.transpose() * b * tmp.transpose(); - } - else - AminusSigmaB = AminusSigmaB - sigma * - b.transpose() * b; - } - else - AminusSigmaB = AminusSigmaB - sigma * b; - } - else - { - SparseMatrix sigmat (n, n, n); - - // Create sigma * speye(n,n) - sigmat.xcidx (0) = 0; - for (octave_idx_type i = 0; i < n; i++) - { - sigmat.xdata(i) = sigma; - sigmat.xridx(i) = i; - sigmat.xcidx(i+1) = i + 1; - } - - AminusSigmaB = AminusSigmaB - sigmat; - } - - SparseLU fact (AminusSigmaB); - - L = fact.L (); - U = fact.U (); - const octave_idx_type *P2 = fact.row_perm (); - const octave_idx_type *Q2 = fact.col_perm (); - - for (octave_idx_type j = 0; j < n; j++) - { - P[j] = P2[j]; - Q[j] = Q2[j]; - } - - // Test condition number of LU decomposition - double minU = octave_NaN; - double maxU = octave_NaN; - for (octave_idx_type j = 0; j < n; j++) - { - double d = 0.; - if (U.xcidx(j+1) > U.xcidx(j) && - U.xridx (U.xcidx(j+1)-1) == j) - d = std::abs (U.xdata (U.xcidx(j+1)-1)); - - if (xisnan (minU) || d < minU) - minU = d; - - if (xisnan (maxU) || d > maxU) - maxU = d; - } - - double rcond = (minU / maxU); - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - (*current_liboctave_warning_handler) - ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); - (*current_liboctave_warning_handler) - (" an eigenvalue. Convergence is not guaranteed"); - } - - return true; -} - -static bool -LuAminusSigmaB (const Matrix &m, const Matrix &b, - bool cholB, const ColumnVector& permB, double sigma, - Matrix &L, Matrix &U, octave_idx_type *P, - octave_idx_type *Q) -{ - bool have_b = ! b.is_empty (); - octave_idx_type n = m.cols(); - - // Caclulate LU decomposition of 'A - sigma * B' - Matrix AminusSigmaB (m); - - if (have_b) - { - if (cholB) - { - Matrix tmp = sigma * b.transpose() * b; - const double *pB = permB.fortran_vec(); - double *p = AminusSigmaB.fortran_vec(); - - if (permB.length()) - { - for (octave_idx_type j = 0; - j < b.cols(); j++) - for (octave_idx_type i = 0; - i < b.rows(); i++) - *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]), - static_cast<octave_idx_type>(pB[j])); - } - else - AminusSigmaB = AminusSigmaB - tmp; - } - else - AminusSigmaB = AminusSigmaB - sigma * b; - } - else - { - double *p = AminusSigmaB.fortran_vec(); - - for (octave_idx_type i = 0; i < n; i++) - p[i*(n+1)] -= sigma; - } - - LU fact (AminusSigmaB); - - L = fact.P().transpose() * fact.L (); - U = fact.U (); - for (octave_idx_type j = 0; j < n; j++) - P[j] = Q[j] = j; - - // Test condition number of LU decomposition - double minU = octave_NaN; - double maxU = octave_NaN; - for (octave_idx_type j = 0; j < n; j++) - { - double d = std::abs (U.xelem(j,j)); - if (xisnan (minU) || d < minU) - minU = d; - - if (xisnan (maxU) || d > maxU) - maxU = d; - } - - double rcond = (minU / maxU); - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - (*current_liboctave_warning_handler) - ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); - (*current_liboctave_warning_handler) - (" an eigenvalue. Convergence is not guaranteed"); - } - - return true; -} - -static bool -LuAminusSigmaB (const SparseComplexMatrix &m, const SparseComplexMatrix &b, - bool cholB, const ColumnVector& permB, Complex sigma, - SparseComplexMatrix &L, SparseComplexMatrix &U, - octave_idx_type *P, octave_idx_type *Q) -{ - bool have_b = ! b.is_empty (); - octave_idx_type n = m.rows(); - - // Caclulate LU decomposition of 'A - sigma * B' - SparseComplexMatrix AminusSigmaB (m); - - if (have_b) - { - if (cholB) - { - if (permB.length()) - { - SparseMatrix tmp(n,n,n); - for (octave_idx_type i = 0; i < n; i++) - { - tmp.xcidx(i) = i; - tmp.xridx(i) = - static_cast<octave_idx_type>(permB(i)); - tmp.xdata(i) = 1; - } - tmp.xcidx(n) = n; - - AminusSigmaB = AminusSigmaB - tmp * b.hermitian() * b * - tmp.transpose() * sigma; - } - else - AminusSigmaB = AminusSigmaB - sigma * b.hermitian() * b; - } - else - AminusSigmaB = AminusSigmaB - sigma * b; - } - else - { - SparseComplexMatrix sigmat (n, n, n); - - // Create sigma * speye(n,n) - sigmat.xcidx (0) = 0; - for (octave_idx_type i = 0; i < n; i++) - { - sigmat.xdata(i) = sigma; - sigmat.xridx(i) = i; - sigmat.xcidx(i+1) = i + 1; - } - - AminusSigmaB = AminusSigmaB - sigmat; - } - - SparseComplexLU fact (AminusSigmaB); - - L = fact.L (); - U = fact.U (); - const octave_idx_type *P2 = fact.row_perm (); - const octave_idx_type *Q2 = fact.col_perm (); - - for (octave_idx_type j = 0; j < n; j++) - { - P[j] = P2[j]; - Q[j] = Q2[j]; - } - - // Test condition number of LU decomposition - double minU = octave_NaN; - double maxU = octave_NaN; - for (octave_idx_type j = 0; j < n; j++) - { - double d = 0.; - if (U.xcidx(j+1) > U.xcidx(j) && - U.xridx (U.xcidx(j+1)-1) == j) - d = std::abs (U.xdata (U.xcidx(j+1)-1)); - - if (xisnan (minU) || d < minU) - minU = d; - - if (xisnan (maxU) || d > maxU) - maxU = d; - } - - double rcond = (minU / maxU); - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - (*current_liboctave_warning_handler) - ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); - (*current_liboctave_warning_handler) - (" an eigenvalue. Convergence is not guaranteed"); - } - - return true; -} - -static bool -LuAminusSigmaB (const ComplexMatrix &m, const ComplexMatrix &b, - bool cholB, const ColumnVector& permB, Complex sigma, - ComplexMatrix &L, ComplexMatrix &U, octave_idx_type *P, - octave_idx_type *Q) -{ - bool have_b = ! b.is_empty (); - octave_idx_type n = m.cols(); - - // Caclulate LU decomposition of 'A - sigma * B' - ComplexMatrix AminusSigmaB (m); - - if (have_b) - { - if (cholB) - { - ComplexMatrix tmp = sigma * b.hermitian() * b; - const double *pB = permB.fortran_vec(); - Complex *p = AminusSigmaB.fortran_vec(); - - if (permB.length()) - { - for (octave_idx_type j = 0; - j < b.cols(); j++) - for (octave_idx_type i = 0; - i < b.rows(); i++) - *p++ -= tmp.xelem (static_cast<octave_idx_type>(pB[i]), - static_cast<octave_idx_type>(pB[j])); - } - else - AminusSigmaB = AminusSigmaB - tmp; - } - else - AminusSigmaB = AminusSigmaB - sigma * b; - } - else - { - Complex *p = AminusSigmaB.fortran_vec(); - - for (octave_idx_type i = 0; i < n; i++) - p[i*(n+1)] -= sigma; - } - - ComplexLU fact (AminusSigmaB); - - L = fact.P().transpose() * fact.L (); - U = fact.U (); - for (octave_idx_type j = 0; j < n; j++) - P[j] = Q[j] = j; - - // Test condition number of LU decomposition - double minU = octave_NaN; - double maxU = octave_NaN; - for (octave_idx_type j = 0; j < n; j++) - { - double d = std::abs (U.xelem(j,j)); - if (xisnan (minU) || d < minU) - minU = d; - - if (xisnan (maxU) || d > maxU) - maxU = d; - } - - double rcond = (minU / maxU); - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) - { - (*current_liboctave_warning_handler) - ("eigs: 'A - sigma*B' is singular, indicating sigma is exactly"); - (*current_liboctave_warning_handler) - (" an eigenvalue. Convergence is not guaranteed"); - } - - return true; -} - -template <class M> -octave_idx_type -EigsRealSymmetricMatrix (const M& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const M& _b, - ColumnVector &permB, ColumnVector &resid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) -{ - M b(_b); - octave_idx_type n = m.cols (); - octave_idx_type mode = 1; - bool have_b = ! b.is_empty(); - bool note3 = false; - char bmat = 'I'; - double sigma = 0.; - M bt; - - if (m.rows() != m.cols()) - { - (*current_liboctave_error_handler) ("eigs: A must be square"); - return -1; - } - if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) - { - (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); - return -1; - } - - if (resid.is_empty()) - { - std::string rand_dist = octave_rand::distribution(); - octave_rand::distribution("uniform"); - resid = ColumnVector (octave_rand::vector(n)); - octave_rand::distribution(rand_dist); - } - - if (p < 0) - { - p = k * 2; - - if (p < 20) - p = 20; - - if (p > n - 1) - p = n - 1 ; - } - else if (p <= k || p > n) - { - (*current_liboctave_error_handler) - ("eigs: opts.p must be between k and n"); - return -1; - } - - if (k > n ) - { - (*current_liboctave_error_handler) - ("eigs: Too many eigenvalues to extract (k >= n).\n" - " Use 'eig(full(A))' instead"); - return -1; - } - - if (have_b && cholB && permB.length() != 0) - { - // Check the we really have a permutation vector - if (permB.length() != n) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } - } - - if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && - typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && - typ != "SI") - { - (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); - return -1; - } - - if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR") - { - (*current_liboctave_error_handler) - ("eigs: invalid sigma value for real symmetric problem"); - return -1; - } - - if (have_b) - { - // See Note 3 dsaupd - note3 = true; - if (cholB) - { - bt = b; - b = b.transpose(); - if (permB.length() == 0) - { - permB = ColumnVector(n); - for (octave_idx_type i = 0; i < n; i++) - permB(i) = i; - } - } - else - { - if (! make_cholb(b, bt, permB)) - { - (*current_liboctave_error_handler) - ("eigs: The matrix B is not positive definite"); - return -1; - } - } - } - - Array<octave_idx_type> ip (11); - octave_idx_type *iparam = ip.fortran_vec (); - - ip(0) = 1; //ishift - ip(1) = 0; // ip(1) not referenced - ip(2) = maxit; // mxiter, maximum number of iterations - ip(3) = 1; // NB blocksize in recurrence - ip(4) = 0; // nconv, number of Ritz values that satisfy convergence - ip(5) = 0; //ip(5) not referenced - ip(6) = mode; // mode - ip(7) = 0; - ip(8) = 0; - ip(9) = 0; - ip(10) = 0; - // ip(7) to ip(10) return values - - Array<octave_idx_type> iptr (14); - octave_idx_type *ipntr = iptr.fortran_vec (); - - octave_idx_type ido = 0; - int iter = 0; - octave_idx_type lwork = p * (p + 8); - - OCTAVE_LOCAL_BUFFER (double, v, n * p); - OCTAVE_LOCAL_BUFFER (double, workl, lwork); - OCTAVE_LOCAL_BUFFER (double, workd, 3 * n); - double *presid = resid.fortran_vec (); - - do - { - F77_FUNC (dsaupd, DSAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dsaupd"); - return -1; - } - - if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } - - if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - Matrix mtmp (n,1); - for (octave_idx_type i = 0; i < n; i++) - mtmp(i,0) = workd[i + iptr(0) - 1]; - - mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); - - for (octave_idx_type i = 0; i < n; i++) - workd[i+iptr(1)-1] = mtmp(i,0); - } - else if (!vector_product (m, workd + iptr(0) - 1, - workd + iptr(1) - 1)) - break; - } - else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } - } - while (1); - - octave_idx_type info2; - - // We have a problem in that the size of the C++ bool - // type relative to the fortran logical type. It appears - // that fortran uses 4-bytes per logical and C++ 1-byte - // per bool, though this might be system dependent. As - // long as the HOWMNY arg is not "S", the logical array - // is just workspace for ARPACK, so use int type to - // avoid problems. - Array<int> s (p); - int *sel = s.fortran_vec (); - - eig_vec.resize (n, k); - double *z = eig_vec.fortran_vec (); - - eig_val.resize (k); - double *d = eig_val.fortran_vec (); - - F77_FUNC (dseupd, DSEUPD) - (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, - F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) - F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dseupd"); - return -1; - } - else - { - if (info2 == 0) - { - octave_idx_type k2 = k / 2; - if (typ != "SM" && typ != "BE") - { - for (octave_idx_type i = 0; i < k2; i++) - { - double dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - } - - if (rvec) - { - if (typ != "SM" && typ != "BE") - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - } - - if (note3) - eig_vec = ltsolve(b, permB, eig_vec); - } - } - else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dseupd", info2); - return -1; - } - } - - return ip(4); -} - -template <class M> -octave_idx_type -EigsRealSymmetricMatrixShift (const M& m, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const M& _b, - ColumnVector &permB, ColumnVector &resid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) -{ - M b(_b); - octave_idx_type n = m.cols (); - octave_idx_type mode = 3; - bool have_b = ! b.is_empty(); - std::string typ = "LM"; - - if (m.rows() != m.cols()) - { - (*current_liboctave_error_handler) ("eigs: A must be square"); - return -1; - } - if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) - { - (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); - return -1; - } - - // FIXME: The "SM" type for mode 1 seems unstable though faster!! - //if (! std::abs (sigma)) - // return EigsRealSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val, - // _b, permB, resid, os, tol, rvec, cholB, - // disp, maxit); - - if (resid.is_empty()) - { - std::string rand_dist = octave_rand::distribution(); - octave_rand::distribution("uniform"); - resid = ColumnVector (octave_rand::vector(n)); - octave_rand::distribution(rand_dist); - } - - if (p < 0) - { - p = k * 2; - - if (p < 20) - p = 20; - - if (p > n - 1) - p = n - 1 ; - } - else if (p <= k || p > n) - { - (*current_liboctave_error_handler) - ("eigs: opts.p must be between k and n"); - return -1; - } - - if (k > n ) - { - (*current_liboctave_error_handler) - ("eigs: Too many eigenvalues to extract (k >= n).\n" - " Use 'eig(full(A))' instead"); - return -1; - } - - if (have_b && cholB && permB.length() != 0) - { - // Check the we really have a permutation vector - if (permB.length() != n) - { - (*current_liboctave_error_handler) ("eigs: permB vector invalid"); - return -1; - } - else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } - } - - char bmat = 'I'; - if (have_b) - bmat = 'G'; - - Array<octave_idx_type> ip (11); - octave_idx_type *iparam = ip.fortran_vec (); - - ip(0) = 1; //ishift - ip(1) = 0; // ip(1) not referenced - ip(2) = maxit; // mxiter, maximum number of iterations - ip(3) = 1; // NB blocksize in recurrence - ip(4) = 0; // nconv, number of Ritz values that satisfy convergence - ip(5) = 0; //ip(5) not referenced - ip(6) = mode; // mode - ip(7) = 0; - ip(8) = 0; - ip(9) = 0; - ip(10) = 0; - // ip(7) to ip(10) return values - - Array<octave_idx_type> iptr (14); - octave_idx_type *ipntr = iptr.fortran_vec (); - - octave_idx_type ido = 0; - int iter = 0; - M L, U; - - OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows())); - OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols())); - - if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q)) - return -1; - - octave_idx_type lwork = p * (p + 8); - - OCTAVE_LOCAL_BUFFER (double, v, n * p); - OCTAVE_LOCAL_BUFFER (double, workl, lwork); - OCTAVE_LOCAL_BUFFER (double, workd, 3 * n); - double *presid = resid.fortran_vec (); - - do - { - F77_FUNC (dsaupd, DSAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dsaupd"); - return -1; - } - - if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } - - if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - if (ido == -1) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - vector_product (m, workd+iptr(0)-1, dtmp); - - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = dtmp[P[i]]; - - lusolve (L, U, tmp); - - double *ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - else if (ido == 2) - vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); - else - { - double *ip2 = workd+iptr(2)-1; - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - else - { - if (ido == 2) - { - for (octave_idx_type i = 0; i < n; i++) - workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1]; - } - else - { - double *ip2 = workd+iptr(0)-1; - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - } - else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } - } - while (1); - - octave_idx_type info2; - - // We have a problem in that the size of the C++ bool - // type relative to the fortran logical type. It appears - // that fortran uses 4-bytes per logical and C++ 1-byte - // per bool, though this might be system dependent. As - // long as the HOWMNY arg is not "S", the logical array - // is just workspace for ARPACK, so use int type to - // avoid problems. - Array<int> s (p); - int *sel = s.fortran_vec (); - - eig_vec.resize (n, k); - double *z = eig_vec.fortran_vec (); - - eig_val.resize (k); - double *d = eig_val.fortran_vec (); - - F77_FUNC (dseupd, DSEUPD) - (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, - F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), - k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2 - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dseupd"); - return -1; - } - else - { - if (info2 == 0) - { - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - double dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - } - } - else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dseupd", info2); - return -1; - } - } - - return ip(4); -} - -octave_idx_type -EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n, - const std::string &_typ, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, ColumnVector &resid, - std::ostream& os, double tol, int rvec, bool cholB, - int disp, int maxit) -{ - std::string typ (_typ); - bool have_sigma = (sigma ? true : false); - char bmat = 'I'; - octave_idx_type mode = 1; - int err = 0; - - if (resid.is_empty()) - { - std::string rand_dist = octave_rand::distribution(); - octave_rand::distribution("uniform"); - resid = ColumnVector (octave_rand::vector(n)); - octave_rand::distribution(rand_dist); - } - - if (p < 0) - { - p = k * 2; - - if (p < 20) - p = 20; - - if (p > n - 1) - p = n - 1 ; - } - else if (p <= k || p > n) - { - (*current_liboctave_error_handler) - ("eigs: opts.p must be between k and n"); - return -1; - } - - if (k > n ) - { - (*current_liboctave_error_handler) - ("eigs: Too many eigenvalues to extract (k >= n).\n" - " Use 'eig(full(A))' instead"); - return -1; - } - - if (! have_sigma) - { - if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && - typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && - typ != "SI") - (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); - - if (typ == "LI" || typ == "SI" || typ == "LR" || typ == "SR") - { - (*current_liboctave_error_handler) - ("eigs: invalid sigma value for real symmetric problem"); - return -1; - } - - if (typ == "SM") - { - typ = "LM"; - sigma = 0.; - mode = 3; - } - } - else if (! std::abs (sigma)) - typ = "SM"; - else - { - typ = "LM"; - mode = 3; - } - - Array<octave_idx_type> ip (11); - octave_idx_type *iparam = ip.fortran_vec (); - - ip(0) = 1; //ishift - ip(1) = 0; // ip(1) not referenced - ip(2) = maxit; // mxiter, maximum number of iterations - ip(3) = 1; // NB blocksize in recurrence - ip(4) = 0; // nconv, number of Ritz values that satisfy convergence - ip(5) = 0; //ip(5) not referenced - ip(6) = mode; // mode - ip(7) = 0; - ip(8) = 0; - ip(9) = 0; - ip(10) = 0; - // ip(7) to ip(10) return values - - Array<octave_idx_type> iptr (14); - octave_idx_type *ipntr = iptr.fortran_vec (); - - octave_idx_type ido = 0; - int iter = 0; - octave_idx_type lwork = p * (p + 8); - - OCTAVE_LOCAL_BUFFER (double, v, n * p); - OCTAVE_LOCAL_BUFFER (double, workl, lwork); - OCTAVE_LOCAL_BUFFER (double, workd, 3 * n); - double *presid = resid.fortran_vec (); - - do - { - F77_FUNC (dsaupd, DSAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dsaupd"); - return -1; - } - - if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } - - - if (ido == -1 || ido == 1 || ido == 2) - { - double *ip2 = workd + iptr(0) - 1; - ColumnVector x(n); - - for (octave_idx_type i = 0; i < n; i++) - x(i) = *ip2++; - - ColumnVector y = fun (x, err); - - if (err) - return false; - - ip2 = workd + iptr(1) - 1; - for (octave_idx_type i = 0; i < n; i++) - *ip2++ = y(i); - } - else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } - } - while (1); - - octave_idx_type info2; - - // We have a problem in that the size of the C++ bool - // type relative to the fortran logical type. It appears - // that fortran uses 4-bytes per logical and C++ 1-byte - // per bool, though this might be system dependent. As - // long as the HOWMNY arg is not "S", the logical array - // is just workspace for ARPACK, so use int type to - // avoid problems. - Array<int> s (p); - int *sel = s.fortran_vec (); - - eig_vec.resize (n, k); - double *z = eig_vec.fortran_vec (); - - eig_val.resize (k); - double *d = eig_val.fortran_vec (); - - F77_FUNC (dseupd, DSEUPD) - (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, - F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), - k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, info2 - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dseupd"); - return -1; - } - else - { - if (info2 == 0) - { - octave_idx_type k2 = k / 2; - if (typ != "SM" && typ != "BE") - { - for (octave_idx_type i = 0; i < k2; i++) - { - double dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - } - - if (rvec) - { - if (typ != "SM" && typ != "BE") - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - } - } - } - else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dseupd", info2); - return -1; - } - } - - return ip(4); -} - -template <class M> -octave_idx_type -EigsRealNonSymmetricMatrix (const M& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const M& _b, - ColumnVector &permB, ColumnVector &resid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) -{ - M b(_b); - octave_idx_type n = m.cols (); - octave_idx_type mode = 1; - bool have_b = ! b.is_empty(); - bool note3 = false; - char bmat = 'I'; - double sigmar = 0.; - double sigmai = 0.; - M bt; - - if (m.rows() != m.cols()) - { - (*current_liboctave_error_handler) ("eigs: A must be square"); - return -1; - } - if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) - { - (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); - return -1; - } - - if (resid.is_empty()) - { - std::string rand_dist = octave_rand::distribution(); - octave_rand::distribution("uniform"); - resid = ColumnVector (octave_rand::vector(n)); - octave_rand::distribution(rand_dist); - } - - if (p < 0) - { - p = k * 2 + 1; - - if (p < 20) - p = 20; - - if (p > n - 1) - p = n - 1 ; - } - else if (p < k || p > n) - { - (*current_liboctave_error_handler) - ("eigs: opts.p must be between k+1 and n"); - return -1; - } - - if (k > n - 1) - { - (*current_liboctave_error_handler) - ("eigs: Too many eigenvalues to extract (k >= n-1).\n" - " Use 'eig(full(A))' instead"); - return -1; - } - - if (have_b && cholB && permB.length() != 0) - { - // Check the we really have a permutation vector - if (permB.length() != n) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } - } - - if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && - typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && - typ != "SI") - { - (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); - return -1; - } - - if (typ == "LA" || typ == "SA" || typ == "BE") - { - (*current_liboctave_error_handler) - ("eigs: invalid sigma value for unsymmetric problem"); - return -1; - } - - if (have_b) - { - // See Note 3 dsaupd - note3 = true; - if (cholB) - { - bt = b; - b = b.transpose(); - if (permB.length() == 0) - { - permB = ColumnVector(n); - for (octave_idx_type i = 0; i < n; i++) - permB(i) = i; - } - } - else - { - if (! make_cholb(b, bt, permB)) - { - (*current_liboctave_error_handler) - ("eigs: The matrix B is not positive definite"); - return -1; - } - } - } - - Array<octave_idx_type> ip (11); - octave_idx_type *iparam = ip.fortran_vec (); - - ip(0) = 1; //ishift - ip(1) = 0; // ip(1) not referenced - ip(2) = maxit; // mxiter, maximum number of iterations - ip(3) = 1; // NB blocksize in recurrence - ip(4) = 0; // nconv, number of Ritz values that satisfy convergence - ip(5) = 0; //ip(5) not referenced - ip(6) = mode; // mode - ip(7) = 0; - ip(8) = 0; - ip(9) = 0; - ip(10) = 0; - // ip(7) to ip(10) return values - - Array<octave_idx_type> iptr (14); - octave_idx_type *ipntr = iptr.fortran_vec (); - - octave_idx_type ido = 0; - int iter = 0; - octave_idx_type lwork = 3 * p * (p + 2); - - OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1)); - OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1); - OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1); - double *presid = resid.fortran_vec (); - - do - { - F77_FUNC (dnaupd, DNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dnaupd"); - return -1; - } - - if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } - - if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - Matrix mtmp (n,1); - for (octave_idx_type i = 0; i < n; i++) - mtmp(i,0) = workd[i + iptr(0) - 1]; - - mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); - - for (octave_idx_type i = 0; i < n; i++) - workd[i+iptr(1)-1] = mtmp(i,0); - } - else if (!vector_product (m, workd + iptr(0) - 1, - workd + iptr(1) - 1)) - break; - } - else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dnaupd", info); - return -1; - } - break; - } - } - while (1); - - octave_idx_type info2; - - // We have a problem in that the size of the C++ bool - // type relative to the fortran logical type. It appears - // that fortran uses 4-bytes per logical and C++ 1-byte - // per bool, though this might be system dependent. As - // long as the HOWMNY arg is not "S", the logical array - // is just workspace for ARPACK, so use int type to - // avoid problems. - Array<int> s (p); - int *sel = s.fortran_vec (); - - Matrix eig_vec2 (n, k + 1); - double *z = eig_vec2.fortran_vec (); - - OCTAVE_LOCAL_BUFFER (double, dr, k + 1); - OCTAVE_LOCAL_BUFFER (double, di, k + 1); - OCTAVE_LOCAL_BUFFER (double, workev, 3 * p); - - F77_FUNC (dneupd, DNEUPD) - (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, - sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) - F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dneupd"); - return -1; - } - else - { - eig_val.resize (k+1); - Complex *d = eig_val.fortran_vec (); - - if (info2 == 0) - { - octave_idx_type jj = 0; - for (octave_idx_type i = 0; i < k+1; i++) - { - if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) - jj++; - else - d [i-jj] = Complex (dr[i], di[i]); - } - if (jj == 0 && !rvec) - for (octave_idx_type i = 0; i < k; i++) - d[i] = d[i+1]; - - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - Complex dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - eig_val.resize(k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - - eig_vec.resize (n, k); - octave_idx_type i = 0; - while (i < k) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (i+1) * n; - if (std::imag(eig_val(i)) == 0) - { - for (octave_idx_type j = 0; j < n; j++) - eig_vec(j,i) = - Complex(z[j+off1],0.); - i++; - } - else - { - for (octave_idx_type j = 0; j < n; j++) - { - eig_vec(j,i) = - Complex(z[j+off1],z[j+off2]); - if (i < k - 1) - eig_vec(j,i+1) = - Complex(z[j+off1],-z[j+off2]); - } - i+=2; - } - } - - if (note3) - eig_vec = ltsolve(M (b), permB, eig_vec); - } - } - else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dneupd", info2); - return -1; - } - } - - return ip(4); -} - -template <class M> -octave_idx_type -EigsRealNonSymmetricMatrixShift (const M& m, double sigmar, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const M& _b, - ColumnVector &permB, ColumnVector &resid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) -{ - M b(_b); - octave_idx_type n = m.cols (); - octave_idx_type mode = 3; - bool have_b = ! b.is_empty(); - std::string typ = "LM"; - double sigmai = 0.; - - if (m.rows() != m.cols()) - { - (*current_liboctave_error_handler) ("eigs: A must be square"); - return -1; - } - if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) - { - (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); - return -1; - } - - // FIXME: The "SM" type for mode 1 seems unstable though faster!! - //if (! std::abs (sigmar)) - // return EigsRealNonSymmetricMatrix (m, "SM", k, p, info, eig_vec, eig_val, - // _b, permB, resid, os, tol, rvec, cholB, - // disp, maxit); - - if (resid.is_empty()) - { - std::string rand_dist = octave_rand::distribution(); - octave_rand::distribution("uniform"); - resid = ColumnVector (octave_rand::vector(n)); - octave_rand::distribution(rand_dist); - } - - if (p < 0) - { - p = k * 2 + 1; - - if (p < 20) - p = 20; - - if (p > n - 1) - p = n - 1 ; - } - else if (p < k || p > n) - { - (*current_liboctave_error_handler) - ("eigs: opts.p must be between k+1 and n"); - return -1; - } - - if (k > n - 1) - { - (*current_liboctave_error_handler) - ("eigs: Too many eigenvalues to extract (k >= n-1).\n" - " Use 'eig(full(A))' instead"); - return -1; - } - - if (have_b && cholB && permB.length() != 0) - { - // Check that we really have a permutation vector - if (permB.length() != n) - { - (*current_liboctave_error_handler) ("eigs: permB vector invalid"); - return -1; - } - else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } - } - - char bmat = 'I'; - if (have_b) - bmat = 'G'; - - Array<octave_idx_type> ip (11); - octave_idx_type *iparam = ip.fortran_vec (); - - ip(0) = 1; //ishift - ip(1) = 0; // ip(1) not referenced - ip(2) = maxit; // mxiter, maximum number of iterations - ip(3) = 1; // NB blocksize in recurrence - ip(4) = 0; // nconv, number of Ritz values that satisfy convergence - ip(5) = 0; //ip(5) not referenced - ip(6) = mode; // mode - ip(7) = 0; - ip(8) = 0; - ip(9) = 0; - ip(10) = 0; - // ip(7) to ip(10) return values - - Array<octave_idx_type> iptr (14); - octave_idx_type *ipntr = iptr.fortran_vec (); - - octave_idx_type ido = 0; - int iter = 0; - M L, U; - - OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows())); - OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols())); - - if (! LuAminusSigmaB(m, b, cholB, permB, sigmar, L, U, P, Q)) - return -1; - - octave_idx_type lwork = 3 * p * (p + 2); - - OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1)); - OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1); - OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1); - double *presid = resid.fortran_vec (); - - do - { - F77_FUNC (dnaupd, DNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dsaupd"); - return -1; - } - - if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } - - if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - if (ido == -1) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - vector_product (m, workd+iptr(0)-1, dtmp); - - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = dtmp[P[i]]; - - lusolve (L, U, tmp); - - double *ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - else if (ido == 2) - vector_product (b, workd+iptr(0)-1, workd+iptr(1)-1); - else - { - double *ip2 = workd+iptr(2)-1; - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - else - { - if (ido == 2) - { - for (octave_idx_type i = 0; i < n; i++) - workd[iptr(0) + i - 1] = workd[iptr(1) + i - 1]; - } - else - { - double *ip2 = workd+iptr(0)-1; - Matrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - } - else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } - } - while (1); - - octave_idx_type info2; - - // We have a problem in that the size of the C++ bool - // type relative to the fortran logical type. It appears - // that fortran uses 4-bytes per logical and C++ 1-byte - // per bool, though this might be system dependent. As - // long as the HOWMNY arg is not "S", the logical array - // is just workspace for ARPACK, so use int type to - // avoid problems. - Array<int> s (p); - int *sel = s.fortran_vec (); - - Matrix eig_vec2 (n, k + 1); - double *z = eig_vec2.fortran_vec (); - - OCTAVE_LOCAL_BUFFER (double, dr, k + 1); - OCTAVE_LOCAL_BUFFER (double, di, k + 1); - OCTAVE_LOCAL_BUFFER (double, workev, 3 * p); - - F77_FUNC (dneupd, DNEUPD) - (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, - sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) - F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dneupd"); - return -1; - } - else - { - eig_val.resize (k+1); - Complex *d = eig_val.fortran_vec (); - - if (info2 == 0) - { - octave_idx_type jj = 0; - for (octave_idx_type i = 0; i < k+1; i++) - { - if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) - jj++; - else - d [i-jj] = Complex (dr[i], di[i]); - } - if (jj == 0 && !rvec) - for (octave_idx_type i = 0; i < k; i++) - d[i] = d[i+1]; - - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - Complex dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - eig_val.resize(k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - - eig_vec.resize (n, k); - octave_idx_type i = 0; - while (i < k) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (i+1) * n; - if (std::imag(eig_val(i)) == 0) - { - for (octave_idx_type j = 0; j < n; j++) - eig_vec(j,i) = - Complex(z[j+off1],0.); - i++; - } - else - { - for (octave_idx_type j = 0; j < n; j++) - { - eig_vec(j,i) = - Complex(z[j+off1],z[j+off2]); - if (i < k - 1) - eig_vec(j,i+1) = - Complex(z[j+off1],-z[j+off2]); - } - i+=2; - } - } - } - } - else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dneupd", info2); - return -1; - } - } - - return ip(4); -} - -octave_idx_type -EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n, - const std::string &_typ, double sigmar, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, ColumnVector &resid, - std::ostream& os, double tol, int rvec, bool cholB, - int disp, int maxit) -{ - std::string typ (_typ); - bool have_sigma = (sigmar ? true : false); - char bmat = 'I'; - double sigmai = 0.; - octave_idx_type mode = 1; - int err = 0; - - if (resid.is_empty()) - { - std::string rand_dist = octave_rand::distribution(); - octave_rand::distribution("uniform"); - resid = ColumnVector (octave_rand::vector(n)); - octave_rand::distribution(rand_dist); - } - - if (p < 0) - { - p = k * 2 + 1; - - if (p < 20) - p = 20; - - if (p > n - 1) - p = n - 1 ; - } - else if (p < k || p > n) - { - (*current_liboctave_error_handler) - ("eigs: opts.p must be between k+1 and n"); - return -1; - } - - if (k > n - 1) - { - (*current_liboctave_error_handler) - ("eigs: Too many eigenvalues to extract (k >= n-1).\n" - " Use 'eig(full(A))' instead"); - return -1; - } - - - if (! have_sigma) - { - if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && - typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && - typ != "SI") - (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); - - if (typ == "LA" || typ == "SA" || typ == "BE") - { - (*current_liboctave_error_handler) - ("eigs: invalid sigma value for unsymmetric problem"); - return -1; - } - - if (typ == "SM") - { - typ = "LM"; - sigmar = 0.; - mode = 3; - } - } - else if (! std::abs (sigmar)) - typ = "SM"; - else - { - typ = "LM"; - mode = 3; - } - - Array<octave_idx_type> ip (11); - octave_idx_type *iparam = ip.fortran_vec (); - - ip(0) = 1; //ishift - ip(1) = 0; // ip(1) not referenced - ip(2) = maxit; // mxiter, maximum number of iterations - ip(3) = 1; // NB blocksize in recurrence - ip(4) = 0; // nconv, number of Ritz values that satisfy convergence - ip(5) = 0; //ip(5) not referenced - ip(6) = mode; // mode - ip(7) = 0; - ip(8) = 0; - ip(9) = 0; - ip(10) = 0; - // ip(7) to ip(10) return values - - Array<octave_idx_type> iptr (14); - octave_idx_type *ipntr = iptr.fortran_vec (); - - octave_idx_type ido = 0; - int iter = 0; - octave_idx_type lwork = 3 * p * (p + 2); - - OCTAVE_LOCAL_BUFFER (double, v, n * (p + 1)); - OCTAVE_LOCAL_BUFFER (double, workl, lwork + 1); - OCTAVE_LOCAL_BUFFER (double, workd, 3 * n + 1); - double *presid = resid.fortran_vec (); - - do - { - F77_FUNC (dnaupd, DNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dnaupd"); - return -1; - } - - if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } - - if (ido == -1 || ido == 1 || ido == 2) - { - double *ip2 = workd + iptr(0) - 1; - ColumnVector x(n); - - for (octave_idx_type i = 0; i < n; i++) - x(i) = *ip2++; - - ColumnVector y = fun (x, err); - - if (err) - return false; - - ip2 = workd + iptr(1) - 1; - for (octave_idx_type i = 0; i < n; i++) - *ip2++ = y(i); - } - else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } - } - while (1); - - octave_idx_type info2; - - // We have a problem in that the size of the C++ bool - // type relative to the fortran logical type. It appears - // that fortran uses 4-bytes per logical and C++ 1-byte - // per bool, though this might be system dependent. As - // long as the HOWMNY arg is not "S", the logical array - // is just workspace for ARPACK, so use int type to - // avoid problems. - Array<int> s (p); - int *sel = s.fortran_vec (); - - Matrix eig_vec2 (n, k + 1); - double *z = eig_vec2.fortran_vec (); - - OCTAVE_LOCAL_BUFFER (double, dr, k + 1); - OCTAVE_LOCAL_BUFFER (double, di, k + 1); - OCTAVE_LOCAL_BUFFER (double, workev, 3 * p); - - F77_FUNC (dneupd, DNEUPD) - (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar, - sigmai, workev, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) - F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in dneupd"); - return -1; - } - else - { - eig_val.resize (k+1); - Complex *d = eig_val.fortran_vec (); - - if (info2 == 0) - { - octave_idx_type jj = 0; - for (octave_idx_type i = 0; i < k+1; i++) - { - if (dr[i] == 0.0 && di[i] == 0.0 && jj == 0) - jj++; - else - d [i-jj] = Complex (dr[i], di[i]); - } - if (jj == 0 && !rvec) - for (octave_idx_type i = 0; i < k; i++) - d[i] = d[i+1]; - - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - Complex dtmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = dtmp; - } - eig_val.resize(k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (double, dtmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - dtmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = dtmp[j]; - } - - eig_vec.resize (n, k); - octave_idx_type i = 0; - while (i < k) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (i+1) * n; - if (std::imag(eig_val(i)) == 0) - { - for (octave_idx_type j = 0; j < n; j++) - eig_vec(j,i) = - Complex(z[j+off1],0.); - i++; - } - else - { - for (octave_idx_type j = 0; j < n; j++) - { - eig_vec(j,i) = - Complex(z[j+off1],z[j+off2]); - if (i < k - 1) - eig_vec(j,i+1) = - Complex(z[j+off1],-z[j+off2]); - } - i+=2; - } - } - } - } - else - { - (*current_liboctave_error_handler) - ("eigs: error %d in dneupd", info2); - return -1; - } - } - - return ip(4); -} - -template <class M> -octave_idx_type -EigsComplexNonSymmetricMatrix (const M& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const M& _b, - ColumnVector &permB, - ComplexColumnVector &cresid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) -{ - M b(_b); - octave_idx_type n = m.cols (); - octave_idx_type mode = 1; - bool have_b = ! b.is_empty(); - bool note3 = false; - char bmat = 'I'; - Complex sigma = 0.; - M bt; - - if (m.rows() != m.cols()) - { - (*current_liboctave_error_handler) ("eigs: A must be square"); - return -1; - } - if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) - { - (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); - return -1; - } - - if (cresid.is_empty()) - { - std::string rand_dist = octave_rand::distribution(); - octave_rand::distribution("uniform"); - Array<double> rr (octave_rand::vector(n)); - Array<double> ri (octave_rand::vector(n)); - cresid = ComplexColumnVector (n); - for (octave_idx_type i = 0; i < n; i++) - cresid(i) = Complex(rr(i),ri(i)); - octave_rand::distribution(rand_dist); - } - - if (p < 0) - { - p = k * 2 + 1; - - if (p < 20) - p = 20; - - if (p > n - 1) - p = n - 1 ; - } - else if (p < k || p > n) - { - (*current_liboctave_error_handler) - ("eigs: opts.p must be between k+1 and n"); - return -1; - } - - if (k > n - 1) - { - (*current_liboctave_error_handler) - ("eigs: Too many eigenvalues to extract (k >= n-1).\n" - " Use 'eig(full(A))' instead"); - return -1; - } - - if (have_b && cholB && permB.length() != 0) - { - // Check the we really have a permutation vector - if (permB.length() != n) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } - } - - if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && - typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && - typ != "SI") - { - (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); - return -1; - } - - if (typ == "LA" || typ == "SA" || typ == "BE") - { - (*current_liboctave_error_handler) - ("eigs: invalid sigma value for complex problem"); - return -1; - } - - if (have_b) - { - // See Note 3 dsaupd - note3 = true; - if (cholB) - { - bt = b; - b = b.hermitian(); - if (permB.length() == 0) - { - permB = ColumnVector(n); - for (octave_idx_type i = 0; i < n; i++) - permB(i) = i; - } - } - else - { - if (! make_cholb(b, bt, permB)) - { - (*current_liboctave_error_handler) - ("eigs: The matrix B is not positive definite"); - return -1; - } - } - } - - Array<octave_idx_type> ip (11); - octave_idx_type *iparam = ip.fortran_vec (); - - ip(0) = 1; //ishift - ip(1) = 0; // ip(1) not referenced - ip(2) = maxit; // mxiter, maximum number of iterations - ip(3) = 1; // NB blocksize in recurrence - ip(4) = 0; // nconv, number of Ritz values that satisfy convergence - ip(5) = 0; //ip(5) not referenced - ip(6) = mode; // mode - ip(7) = 0; - ip(8) = 0; - ip(9) = 0; - ip(10) = 0; - // ip(7) to ip(10) return values - - Array<octave_idx_type> iptr (14); - octave_idx_type *ipntr = iptr.fortran_vec (); - - octave_idx_type ido = 0; - int iter = 0; - octave_idx_type lwork = p * (3 * p + 5); - - OCTAVE_LOCAL_BUFFER (Complex, v, n * p); - OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); - OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); - OCTAVE_LOCAL_BUFFER (double, rwork, p); - Complex *presid = cresid.fortran_vec (); - - do - { - F77_FUNC (znaupd, ZNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, rwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in znaupd"); - return -1; - } - - if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } - - if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - ComplexMatrix mtmp (n,1); - for (octave_idx_type i = 0; i < n; i++) - mtmp(i,0) = workd[i + iptr(0) - 1]; - mtmp = utsolve(bt, permB, m * ltsolve(b, permB, mtmp)); - for (octave_idx_type i = 0; i < n; i++) - workd[i+iptr(1)-1] = mtmp(i,0); - - } - else if (!vector_product (m, workd + iptr(0) - 1, - workd + iptr(1) - 1)) - break; - } - else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in znaupd", info); - return -1; - } - break; - } - } - while (1); - - octave_idx_type info2; - - // We have a problem in that the size of the C++ bool - // type relative to the fortran logical type. It appears - // that fortran uses 4-bytes per logical and C++ 1-byte - // per bool, though this might be system dependent. As - // long as the HOWMNY arg is not "S", the logical array - // is just workspace for ARPACK, so use int type to - // avoid problems. - Array<int> s (p); - int *sel = s.fortran_vec (); - - eig_vec.resize (n, k); - Complex *z = eig_vec.fortran_vec (); - - eig_val.resize (k+1); - Complex *d = eig_val.fortran_vec (); - - OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p); - - F77_FUNC (zneupd, ZNEUPD) - (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev, - F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), - k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2 - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in zneupd"); - return -1; - } - - if (info2 == 0) - { - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - Complex ctmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = ctmp; - } - eig_val.resize(k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - ctmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = ctmp[j]; - } - - if (note3) - eig_vec = ltsolve(b, permB, eig_vec); - } - } - else - { - (*current_liboctave_error_handler) - ("eigs: error %d in zneupd", info2); - return -1; - } - - return ip(4); -} - -template <class M> -octave_idx_type -EigsComplexNonSymmetricMatrixShift (const M& m, Complex sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const M& _b, - ColumnVector &permB, - ComplexColumnVector &cresid, - std::ostream& os, double tol, int rvec, - bool cholB, int disp, int maxit) -{ - M b(_b); - octave_idx_type n = m.cols (); - octave_idx_type mode = 3; - bool have_b = ! b.is_empty(); - std::string typ = "LM"; - - if (m.rows() != m.cols()) - { - (*current_liboctave_error_handler) ("eigs: A must be square"); - return -1; - } - if (have_b && (m.rows() != b.rows() || m.rows() != b.cols())) - { - (*current_liboctave_error_handler) - ("eigs: B must be square and the same size as A"); - return -1; - } - - // FIXME: The "SM" type for mode 1 seems unstable though faster!! - //if (! std::abs (sigma)) - // return EigsComplexNonSymmetricMatrix (m, "SM", k, p, info, eig_vec, - // eig_val, _b, permB, cresid, os, tol, - // rvec, cholB, disp, maxit); - - if (cresid.is_empty()) - { - std::string rand_dist = octave_rand::distribution(); - octave_rand::distribution("uniform"); - Array<double> rr (octave_rand::vector(n)); - Array<double> ri (octave_rand::vector(n)); - cresid = ComplexColumnVector (n); - for (octave_idx_type i = 0; i < n; i++) - cresid(i) = Complex(rr(i),ri(i)); - octave_rand::distribution(rand_dist); - } - - if (p < 0) - { - p = k * 2 + 1; - - if (p < 20) - p = 20; - - if (p > n - 1) - p = n - 1 ; - } - else if (p < k || p > n) - { - (*current_liboctave_error_handler) - ("eigs: opts.p must be between k+1 and n"); - return -1; - } - - if (k > n - 1) - { - (*current_liboctave_error_handler) - ("eigs: Too many eigenvalues to extract (k >= n-1).\n" - " Use 'eig(full(A))' instead"); - return -1; - } - - if (have_b && cholB && permB.length() != 0) - { - // Check that we really have a permutation vector - if (permB.length() != n) - { - (*current_liboctave_error_handler) ("eigs: permB vector invalid"); - return -1; - } - else - { - Array<bool> checked(n,false); - for (octave_idx_type i = 0; i < n; i++) - { - octave_idx_type bidx = - static_cast<octave_idx_type> (permB(i)); - if (checked(bidx) || bidx < 0 || - bidx >= n || D_NINT (bidx) != bidx) - { - (*current_liboctave_error_handler) - ("eigs: permB vector invalid"); - return -1; - } - } - } - } - - char bmat = 'I'; - if (have_b) - bmat = 'G'; - - Array<octave_idx_type> ip (11); - octave_idx_type *iparam = ip.fortran_vec (); - - ip(0) = 1; //ishift - ip(1) = 0; // ip(1) not referenced - ip(2) = maxit; // mxiter, maximum number of iterations - ip(3) = 1; // NB blocksize in recurrence - ip(4) = 0; // nconv, number of Ritz values that satisfy convergence - ip(5) = 0; //ip(5) not referenced - ip(6) = mode; // mode - ip(7) = 0; - ip(8) = 0; - ip(9) = 0; - ip(10) = 0; - // ip(7) to ip(10) return values - - Array<octave_idx_type> iptr (14); - octave_idx_type *ipntr = iptr.fortran_vec (); - - octave_idx_type ido = 0; - int iter = 0; - M L, U; - - OCTAVE_LOCAL_BUFFER (octave_idx_type, P, (have_b ? b.rows() : m.rows())); - OCTAVE_LOCAL_BUFFER (octave_idx_type, Q, (have_b ? b.cols() : m.cols())); - - if (! LuAminusSigmaB(m, b, cholB, permB, sigma, L, U, P, Q)) - return -1; - - octave_idx_type lwork = p * (3 * p + 5); - - OCTAVE_LOCAL_BUFFER (Complex, v, n * p); - OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); - OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); - OCTAVE_LOCAL_BUFFER (double, rwork, p); - Complex *presid = cresid.fortran_vec (); - - do - { - F77_FUNC (znaupd, ZNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, rwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in znaupd"); - return -1; - } - - if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } - - if (ido == -1 || ido == 1 || ido == 2) - { - if (have_b) - { - if (ido == -1) - { - OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); - - vector_product (m, workd+iptr(0)-1, ctmp); - - ComplexMatrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ctmp[P[i]]; - - lusolve (L, U, tmp); - - Complex *ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - else if (ido == 2) - vector_product (b, workd + iptr(0) - 1, workd + iptr(1) - 1); - else - { - Complex *ip2 = workd+iptr(2)-1; - ComplexMatrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - else - { - if (ido == 2) - { - for (octave_idx_type i = 0; i < n; i++) - workd[iptr(0) + i - 1] = - workd[iptr(1) + i - 1]; - } - else - { - Complex *ip2 = workd+iptr(0)-1; - ComplexMatrix tmp(n, 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp(i,0) = ip2[P[i]]; - - lusolve (L, U, tmp); - - ip2 = workd+iptr(1)-1; - for (octave_idx_type i = 0; i < n; i++) - ip2[Q[i]] = tmp(i,0); - } - } - } - else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } - } - while (1); - - octave_idx_type info2; - - // We have a problem in that the size of the C++ bool - // type relative to the fortran logical type. It appears - // that fortran uses 4-bytes per logical and C++ 1-byte - // per bool, though this might be system dependent. As - // long as the HOWMNY arg is not "S", the logical array - // is just workspace for ARPACK, so use int type to - // avoid problems. - Array<int> s (p); - int *sel = s.fortran_vec (); - - eig_vec.resize (n, k); - Complex *z = eig_vec.fortran_vec (); - - eig_val.resize (k+1); - Complex *d = eig_val.fortran_vec (); - - OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p); - - F77_FUNC (zneupd, ZNEUPD) - (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev, - F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), - k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2 - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in zneupd"); - return -1; - } - - if (info2 == 0) - { - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - Complex ctmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = ctmp; - } - eig_val.resize(k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - ctmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = ctmp[j]; - } - } - } - else - { - (*current_liboctave_error_handler) - ("eigs: error %d in zneupd", info2); - return -1; - } - - return ip(4); -} - -octave_idx_type -EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n, - const std::string &_typ, Complex sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - ComplexColumnVector &cresid, std::ostream& os, - double tol, int rvec, bool cholB, int disp, - int maxit) -{ - std::string typ (_typ); - bool have_sigma = (std::abs(sigma) ? true : false); - char bmat = 'I'; - octave_idx_type mode = 1; - int err = 0; - - if (cresid.is_empty()) - { - std::string rand_dist = octave_rand::distribution(); - octave_rand::distribution("uniform"); - Array<double> rr (octave_rand::vector(n)); - Array<double> ri (octave_rand::vector(n)); - cresid = ComplexColumnVector (n); - for (octave_idx_type i = 0; i < n; i++) - cresid(i) = Complex(rr(i),ri(i)); - octave_rand::distribution(rand_dist); - } - - if (p < 0) - { - p = k * 2 + 1; - - if (p < 20) - p = 20; - - if (p > n - 1) - p = n - 1 ; - } - else if (p < k || p > n) - { - (*current_liboctave_error_handler) - ("eigs: opts.p must be between k+1 and n"); - return -1; - } - - if (k > n - 1) - { - (*current_liboctave_error_handler) - ("eigs: Too many eigenvalues to extract (k >= n-1).\n" - " Use 'eig(full(A))' instead"); - return -1; - } - - if (! have_sigma) - { - if (typ != "LM" && typ != "SM" && typ != "LA" && typ != "SA" && - typ != "BE" && typ != "LR" && typ != "SR" && typ != "LI" && - typ != "SI") - (*current_liboctave_error_handler) - ("eigs: unrecognized sigma value"); - - if (typ == "LA" || typ == "SA" || typ == "BE") - { - (*current_liboctave_error_handler) - ("eigs: invalid sigma value for complex problem"); - return -1; - } - - if (typ == "SM") - { - typ = "LM"; - sigma = 0.; - mode = 3; - } - } - else if (! std::abs (sigma)) - typ = "SM"; - else - { - typ = "LM"; - mode = 3; - } - - Array<octave_idx_type> ip (11); - octave_idx_type *iparam = ip.fortran_vec (); - - ip(0) = 1; //ishift - ip(1) = 0; // ip(1) not referenced - ip(2) = maxit; // mxiter, maximum number of iterations - ip(3) = 1; // NB blocksize in recurrence - ip(4) = 0; // nconv, number of Ritz values that satisfy convergence - ip(5) = 0; //ip(5) not referenced - ip(6) = mode; // mode - ip(7) = 0; - ip(8) = 0; - ip(9) = 0; - ip(10) = 0; - // ip(7) to ip(10) return values - - Array<octave_idx_type> iptr (14); - octave_idx_type *ipntr = iptr.fortran_vec (); - - octave_idx_type ido = 0; - int iter = 0; - octave_idx_type lwork = p * (3 * p + 5); - - OCTAVE_LOCAL_BUFFER (Complex, v, n * p); - OCTAVE_LOCAL_BUFFER (Complex, workl, lwork); - OCTAVE_LOCAL_BUFFER (Complex, workd, 3 * n); - OCTAVE_LOCAL_BUFFER (double, rwork, p); - Complex *presid = cresid.fortran_vec (); - - do - { - F77_FUNC (znaupd, ZNAUPD) - (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str()), 2), - k, tol, presid, p, v, n, iparam, - ipntr, workd, workl, lwork, rwork, info - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in znaupd"); - return -1; - } - - if (disp > 0 && !xisnan(workl[iptr(5)-1])) - { - if (iter++) - { - os << "Iteration " << iter - 1 << - ": a few Ritz values of the " << p << "-by-" << - p << " matrix\n"; - for (int i = 0 ; i < k; i++) - os << " " << workl[iptr(5)+i-1] << "\n"; - } - - // This is a kludge, as ARPACK doesn't give its - // iteration pointer. But as workl[iptr(5)-1] is - // an output value updated at each iteration, setting - // a value in this array to NaN and testing for it - // is a way of obtaining the iteration counter. - if (ido != 99) - workl[iptr(5)-1] = octave_NaN; - } - - if (ido == -1 || ido == 1 || ido == 2) - { - Complex *ip2 = workd + iptr(0) - 1; - ComplexColumnVector x(n); - - for (octave_idx_type i = 0; i < n; i++) - x(i) = *ip2++; - - ComplexColumnVector y = fun (x, err); - - if (err) - return false; - - ip2 = workd + iptr(1) - 1; - for (octave_idx_type i = 0; i < n; i++) - *ip2++ = y(i); - } - else - { - if (info < 0) - { - (*current_liboctave_error_handler) - ("eigs: error %d in dsaupd", info); - return -1; - } - break; - } - } - while (1); - - octave_idx_type info2; - - // We have a problem in that the size of the C++ bool - // type relative to the fortran logical type. It appears - // that fortran uses 4-bytes per logical and C++ 1-byte - // per bool, though this might be system dependent. As - // long as the HOWMNY arg is not "S", the logical array - // is just workspace for ARPACK, so use int type to - // avoid problems. - Array<int> s (p); - int *sel = s.fortran_vec (); - - eig_vec.resize (n, k); - Complex *z = eig_vec.fortran_vec (); - - eig_val.resize (k+1); - Complex *d = eig_val.fortran_vec (); - - OCTAVE_LOCAL_BUFFER (Complex, workev, 2 * p); - - F77_FUNC (zneupd, ZNEUPD) - (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma, workev, - F77_CONST_CHAR_ARG2 (&bmat, 1), n, - F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), - k, tol, presid, p, v, n, iparam, ipntr, workd, workl, lwork, rwork, info2 - F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2)); - - if (f77_exception_encountered) - { - (*current_liboctave_error_handler) - ("eigs: unrecoverable exception encountered in zneupd"); - return -1; - } - - if (info2 == 0) - { - octave_idx_type k2 = k / 2; - for (octave_idx_type i = 0; i < k2; i++) - { - Complex ctmp = d[i]; - d[i] = d[k - i - 1]; - d[k - i - 1] = ctmp; - } - eig_val.resize(k); - - if (rvec) - { - OCTAVE_LOCAL_BUFFER (Complex, ctmp, n); - - for (octave_idx_type i = 0; i < k2; i++) - { - octave_idx_type off1 = i * n; - octave_idx_type off2 = (k - i - 1) * n; - - if (off1 == off2) - continue; - - for (octave_idx_type j = 0; j < n; j++) - ctmp[j] = z[off1 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off1 + j] = z[off2 + j]; - - for (octave_idx_type j = 0; j < n; j++) - z[off2 + j] = ctmp[j]; - } - } - } - else - { - (*current_liboctave_error_handler) - ("eigs: error %d in zneupd", info2); - return -1; - } - - return ip(4); -} - -#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) -extern octave_idx_type -EigsRealSymmetricMatrix (const Matrix& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const Matrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); - -extern octave_idx_type -EigsRealSymmetricMatrix (const SparseMatrix& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const SparseMatrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream& os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); - -extern octave_idx_type -EigsRealSymmetricMatrixShift (const Matrix& m, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const Matrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); - -extern octave_idx_type -EigsRealSymmetricMatrixShift (const SparseMatrix& m, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, Matrix &eig_vec, - ColumnVector &eig_val, const SparseMatrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); - -extern octave_idx_type -EigsRealSymmetricFunc (EigsFunc fun, octave_idx_type n, - const std::string &typ, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - Matrix &eig_vec, ColumnVector &eig_val, - ColumnVector &resid, std::ostream &os, - double tol = DBL_EPSILON, int rvec = 0, - bool cholB = 0, int disp = 0, int maxit = 300); - -extern octave_idx_type -EigsRealNonSymmetricMatrix (const Matrix& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const Matrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); - -extern octave_idx_type -EigsRealNonSymmetricMatrix (const SparseMatrix& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const SparseMatrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); - -extern octave_idx_type -EigsRealNonSymmetricMatrixShift (const Matrix& m, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, const Matrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); - -extern octave_idx_type -EigsRealNonSymmetricMatrixShift (const SparseMatrix& m, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const SparseMatrix& b, - ColumnVector &permB, ColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); - -extern octave_idx_type -EigsRealNonSymmetricFunc (EigsFunc fun, octave_idx_type n, - const std::string &_typ, double sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - ColumnVector &resid, std::ostream& os, - double tol = DBL_EPSILON, int rvec = 0, - bool cholB = 0, int disp = 0, int maxit = 300); - -extern octave_idx_type -EigsComplexNonSymmetricMatrix (const ComplexMatrix& m, const std::string typ, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const ComplexMatrix& b, ColumnVector &permB, - ComplexColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); - -extern octave_idx_type -EigsComplexNonSymmetricMatrix (const SparseComplexMatrix& m, - const std::string typ, octave_idx_type k, - octave_idx_type p, octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const SparseComplexMatrix& b, - ColumnVector &permB, - ComplexColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, int disp = 0, - int maxit = 300); - -extern octave_idx_type -EigsComplexNonSymmetricMatrixShift (const ComplexMatrix& m, Complex sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const ComplexMatrix& b, - ColumnVector &permB, - ComplexColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, - int disp = 0, int maxit = 300); - -extern octave_idx_type -EigsComplexNonSymmetricMatrixShift (const SparseComplexMatrix& m, - Complex sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, - ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - const SparseComplexMatrix& b, - ColumnVector &permB, - ComplexColumnVector &resid, - std::ostream &os, double tol = DBL_EPSILON, - int rvec = 0, bool cholB = 0, - int disp = 0, int maxit = 300); - -extern octave_idx_type -EigsComplexNonSymmetricFunc (EigsComplexFunc fun, octave_idx_type n, - const std::string &_typ, Complex sigma, - octave_idx_type k, octave_idx_type p, - octave_idx_type &info, ComplexMatrix &eig_vec, - ComplexColumnVector &eig_val, - ComplexColumnVector &resid, std::ostream& os, - double tol = DBL_EPSILON, int rvec = 0, - bool cholB = 0, int disp = 0, int maxit = 300); -#endif - -#ifndef _MSC_VER -template static octave_idx_type -lusolve (const SparseMatrix&, const SparseMatrix&, Matrix&); - -template static octave_idx_type -lusolve (const SparseComplexMatrix&, const SparseComplexMatrix&, - ComplexMatrix&); - -template static octave_idx_type -lusolve (const Matrix&, const Matrix&, Matrix&); - -template static octave_idx_type -lusolve (const ComplexMatrix&, const ComplexMatrix&, ComplexMatrix&); - -template static ComplexMatrix -ltsolve (const SparseComplexMatrix&, const ColumnVector&, - const ComplexMatrix&); - -template static Matrix -ltsolve (const SparseMatrix&, const ColumnVector&, const Matrix&); - -template static ComplexMatrix -ltsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&); - -template static Matrix -ltsolve (const Matrix&, const ColumnVector&, const Matrix&); - -template static ComplexMatrix -utsolve (const SparseComplexMatrix&, const ColumnVector&, - const ComplexMatrix&); - -template static Matrix -utsolve (const SparseMatrix&, const ColumnVector&, const Matrix&); - -template static ComplexMatrix -utsolve (const ComplexMatrix&, const ColumnVector&, const ComplexMatrix&); - -template static Matrix -utsolve (const Matrix&, const ColumnVector&, const Matrix&); -#endif - -#endif - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/nonfree/arpack/src/eigs.cc Fri Feb 05 07:21:35 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1463 +0,0 @@ -/* - -Copyright (C) 2005 David Bateman - -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 software 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 software; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "ov.h" -#include "defun-dld.h" -#include "error.h" -#include "gripes.h" -#include "quit.h" -#include "variables.h" -#include "ov-re-sparse.h" -#include "ov-cx-sparse.h" -#include "oct-map.h" -#include "pager.h" -#include "unwind-prot.h" - -#include "eigs-base.cc" - -// Global pointer for user defined function. -static octave_function *eigs_fcn = 0; - -// Have we warned about imaginary values returned from user function? -static bool warned_imaginary = false; - -// Is this a recursive call? -static int call_depth = 0; - -ColumnVector -eigs_func (const ColumnVector &x, int &eigs_error) -{ - ColumnVector retval; - octave_value_list args; - args(0) = x; - - if (eigs_fcn) - { - octave_value_list tmp = eigs_fcn->do_multi_index_op (1, args); - - if (error_state) - { - eigs_error = 1; - gripe_user_supplied_eval ("eigs"); - return retval; - } - - if (tmp.length () && tmp(0).is_defined ()) - { - if (! warned_imaginary && tmp(0).is_complex_type ()) - { - warning ("eigs: ignoring imaginary part returned from user-supplied function"); - warned_imaginary = true; - } - - retval = ColumnVector (tmp(0).vector_value ()); - - if (error_state) - { - eigs_error = 1; - gripe_user_supplied_eval ("eigs"); - } - } - else - { - eigs_error = 1; - gripe_user_supplied_eval ("eigs"); - } - } - - return retval; -} - -ComplexColumnVector -eigs_complex_func (const ComplexColumnVector &x, int &eigs_error) -{ - ComplexColumnVector retval; - octave_value_list args; - args(0) = x; - - if (eigs_fcn) - { - octave_value_list tmp = eigs_fcn->do_multi_index_op (1, args); - - if (error_state) - { - eigs_error = 1; - gripe_user_supplied_eval ("eigs"); - return retval; - } - - if (tmp.length () && tmp(0).is_defined ()) - { - retval = ComplexColumnVector (tmp(0).complex_vector_value ()); - - if (error_state) - { - eigs_error = 1; - gripe_user_supplied_eval ("eigs"); - } - } - else - { - eigs_error = 1; - gripe_user_supplied_eval ("eigs"); - } - } - - return retval; -} - -DEFUN_DLD (eigs, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{d}} = eigs (@var{a})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{k})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{k}, @var{sigma})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{k}, @var{sigma},@var{opts})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b}, @var{k})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b}, @var{k}, @var{sigma})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{a}, @var{b}, @var{k}, @var{sigma}, @var{opts})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{k})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b}, @var{k})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{k}, @var{sigma})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b}, @var{k}, @var{sigma})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{k}, @var{sigma}, @var{opts})\n\ -@deftypefnx {Loadable Function} {@var{d}} = eigs (@var{af}, @var{n}, @var{b}, @var{k}, @var{sigma}, @var{opts})\n\ -@deftypefnx {Loadable Function} {[@var{v}, @var{d}]} = eigs (@var{a}, @dots{})\n\ -@deftypefnx {Loadable Function} {[@var{v}, @var{d}]} = eigs (@var{af}, @var{n}, @dots{})\n\ -@deftypefnx {Loadable Function} {[@var{v}, @var{d}, @var{flag}]} = eigs (@var{a}, @dots{})\n\ -@deftypefnx {Loadable Function} {[@var{v}, @var{d}, @var{flag}]} = eigs (@var{af}, @var{n}, @dots{})\n\ -Calculate a limited number of eigenvalues and eigenvectors of @var{a},\n\ -based on a selection criteria. The number eigenvalues and eigenvectors to\n\ -calculate is given by @var{k} whose default value is 6.\n\ -\n\ -By default @code{eigs} solve the equation\n\ -@iftex\n\ -@tex\n\ -$A \\nu = \\lambda \\nu$\n\ -@end tex\n\ -@end iftex\n\ -@ifinfo\n\ -@code{A * v = lambda * v}\n\ -@end ifinfo\n\ -, where\n\ -@iftex\n\ -@tex\n\ -$\\lambda$ is a scalar representing one of the eigenvalues, and $\\nu$\n\ -@end tex\n\ -@end iftex\n\ -@ifinfo\n\ -@code{lambda} is a scalar representing one of the eigenvalues, and @code{v}\n\ -@end ifinfo\n\ -is the corresponding eigenvector. If given the positive definite matrix\n\ -@var{B} then @code{eigs} solves the general eigenvalue equation\n\ -@iftex\n\ -@tex\n\ -$A \\nu = \\lambda B \\nu$\n\ -@end tex\n\ -@end iftex\n\ -@ifinfo\n\ -@code{A * v = lambda * B * v}\n\ -@end ifinfo\n\ -.\n\ -\n\ -The argument @var{sigma} determines which eigenvalues are returned.\n\ -@var{sigma} can be either a scalar or a string. When @var{sigma} is a scalar,\n\ -the @var{k} eigenvalues closest to @var{sigma} are returned. If @var{sigma}\n\ -is a string, it must have one of the values\n\ -\n\ -@table @asis\n\ -@item 'lm'\n\ -Largest magnitude (default).\n\ -\n\ -@item 'sm'\n\ -Smallest magnitude.\n\ -\n\ -@item 'la'\n\ -Largest Algebraic (valid only for real symmetric problems).\n\ -\n\ -@item 'sa'\n\ -Smallest Algebraic (valid only for real symmetric problems).\n\ -\n\ -@item 'be'\n\ -Both ends, with one more from the high-end if @var{k} is odd (valid only for\n\ -real symmetric problems).\n\ -\n\ -@item 'lr'\n\ -Largest real part (valid only for complex or unsymmetric problems).\n\ -\n\ -@item 'sr'\n\ -Smallest real part (valid only for complex or unsymmetric problems).\n\ -\n\ -@item 'li'\n\ -Largest imaginary part (valid only for complex or unsymmetric problems).\n\ -\n\ -@item 'si'\n\ -Smallest imaginary part (valid only for complex or unsymmetric problems).\n\ -@end table\n\ -\n\ -If @var{opts} is given, it is a structure defining some of the options that\n\ -@code{eigs} should use. The fields of the structure @var{opts} are\n\ -\n\ -@table @code\n\ -@item issym\n\ -If @var{af} is given, then flags whether the function @var{af} defines a\n\ -symmetric problem. It is ignored if @var{a} is given. The default is false.\n\ -\n\ -@item isreal\n\ -If @var{af} is given, then flags whether the function @var{af} defines a\n\ -real problem. It is ignored if @var{a} is given. The default is true.\n\ -\n\ -@item tol\n\ -Defines the required convergence tolerance, given as @code{tol * norm (A)}.\n\ -The default is @code{eps}.\n\ -\n\ -@item maxit\n\ -The maximum number of iterations. The default is 300.\n\ -\n\ -@item p\n\ -The number of Lanzcos basis vectors to use. More vectors will result in\n\ -faster convergence, but a larger amount of memory. The optimal value of 'p'\n\ -is problem dependent and should be in the range @var{k} to @var{n}. The\n\ -default value is @code{2 * @var{k}}.\n\ -\n\ -@item v0\n\ -The starting vector for the computation. The default is to have @sc{Arpack}\n\ -randomly generate a starting vector.\n\ -\n\ -@item disp\n\ -The level of diagnostic printout. If @code{disp} is 0 then there is no\n\ -printout. The default value is 1.\n\ -\n\ -@item cholB\n\ -Flag if @code{chol (@var{b})} is passed rather than @var{b}. The default is\n\ -false.\n\ -\n\ -@item permB\n\ -The permutation vector of the Cholesky factorization of @var{b} if\n\ -@code{cholB} is true. That is @code{chol ( @var{b} (permB, permB))}. The\n\ -default is @code{1:@var{n}}.\n\ -\n\ -@end table\n\ -\n\ -It is also possible to represent @var{a} by a function denoted @var{af}.\n\ -@var{af} must be followed by a scalar argument @var{n} defining the length\n\ -of the vector argument accepted by @var{af}. @var{af} can be passed either\n\ -as an inline function, function handle or as a string. In the case where\n\ -@var{af} is passed as a string, the name of the string defines the function\n\ -to use.\n\ -\n\ -@var{af} is a function of the form @code{function y = af (x), y = @dots{};\n\ -endfunction}, where the required return value of @var{af} is determined by\n\ -the value of @var{sigma}, and are\n\ -\n\ -@table @code\n\ -@item A * x\n\ -If @var{sigma} is not given or is a string other than 'sm'.\n\ -\n\ -@item A \\ x\n\ -If @var{sigma} is 'sm'.\n\ -\n\ -@item (A - sigma * I) \\ x\n\ -for standard eigenvalue problem, where @code{I} is the identity matrix of\n\ -the same size as @code{A}. If @var{sigma} is zero, this reduces the\n\ -@code{A \\ x}.\n\ -\n\ -@item (A - sigma * B) \\ x\n\ -for the general eigenvalue problem.\n\ -@end table\n\ -\n\ -The return arguments of @code{eigs} depends on the number of return\n\ -arguments. With a single return argument, a vector @var{d} of length @var{k}\n\ -is returned, represent the @var{k} eigenvalues that have been found. With two\n\ -return arguments, @var{v} is a @var{n}-by-@var{k} matrix whose columns are\n\ -the @var{k} eigenvectors corresponding to the returned eigenvalues. The\n\ -eigenvalues themselves are then returned in @var{d} in the form of a\n\ -@var{n}-by-@var{k} matrix, where the elements on the diagonal are the\n\ -eigenvalues.\n\ -\n\ -Given a third return argument @var{flag}, @code{eigs} also returns the status\n\ -of the convergence. If @var{flag} is 0, then all eigenvalues have converged,\n\ -otherwise not.\n\ -\n\ -This function is based on the @sc{Arpack} package, written by R Lehoucq,\n\ -K Maschhoff, D Sorensen and C Yang. For more information see\n\ -@url{http://www.caam.rice.edu/software/ARPACK/}.\n\ -\n\ -@end deftypefn\n\ -@seealso{eig, svds}") -{ - octave_value_list retval; -#ifdef HAVE_ARPACK - int nargin = args.length (); - std::string fcn_name; - octave_idx_type n = 0; - octave_idx_type k = 6; - Complex sigma = 0.; - double sigmar, sigmai; - bool have_sigma = false; - std::string typ = "LM"; - Matrix amm, bmm, bmt; - ComplexMatrix acm, bcm, bct; - SparseMatrix asmm, bsmm, bsmt; - SparseComplexMatrix ascm, bscm, bsct; - int b_arg = 0; - bool have_b = false; - bool have_a_fun = false; - bool a_is_complex = false; - bool b_is_complex = false; - bool symmetric = false; - bool cholB = false; - bool a_is_sparse = false; - ColumnVector permB; - int arg_offset = 0; - double tol = DBL_EPSILON; - int maxit = 300; - int disp = 0; - octave_idx_type p = -1; - ColumnVector resid; - ComplexColumnVector cresid; - octave_idx_type info = 1; - char bmat = 'I'; - - warned_imaginary = false; - - unwind_protect::begin_frame ("Feigs"); - - unwind_protect_int (call_depth); - call_depth++; - - if (call_depth > 1) - { - error ("eigs: invalid recursive call"); - if (fcn_name.length()) - clear_function (fcn_name); - unwind_protect::run_frame ("Feigs"); - return retval; - } - - if (nargin == 0) - print_usage (); - else if (args(0).is_function_handle () || args(0).is_inline_function () || - args(0).is_string ()) - { - if (args(0).is_string ()) - { - std::string name = args(0).string_value (); - std::string fname = "function y = "; - fcn_name = unique_symbol_name ("__eigs_fcn_"); - fname.append (fcn_name); - fname.append ("(x) y = "); - eigs_fcn = extract_function (args(0), "eigs", fcn_name, fname, - "; endfunction"); - } - else - eigs_fcn = args(0).function_value (); - - if (!eigs_fcn) - error ("eigs: unknown function"); - - if (nargin < 2) - error ("eigs: incorrect number of arguments"); - else - { - n = args(1).nint_value (); - arg_offset = 1; - have_a_fun = true; - } - } - else - { - if (args(0).is_complex_type ()) - { - if (args(0).is_sparse_type ()) - { - ascm = (args(0).sparse_complex_matrix_value()); - a_is_sparse = true; - } - else - acm = (args(0).complex_matrix_value()); - a_is_complex = true; - symmetric = false; // ARAPACK doesn't special case complex symmetric - } - else - { - if (args(0).is_sparse_type ()) - { - asmm = (args(0).sparse_matrix_value()); - a_is_sparse = true; - symmetric = asmm.is_symmetric(); - } - else - { - amm = (args(0).matrix_value()); - symmetric = amm.is_symmetric(); - } - } - - // Note hold off reading B till later to avoid issues of double - // copies of the matrix if B is full/real while A is complex.. - if (!error_state && nargin > 1 && - !(args(1).is_real_scalar ())) - if (args(1).is_complex_type ()) - { - b_arg = 1+arg_offset; - have_b = true; - bmat = 'G'; - b_is_complex = true; - arg_offset++; - } - else - { - b_arg = 1+arg_offset; - have_b = true; - bmat = 'G'; - arg_offset++; - } - } - - if (!error_state && nargin > (1+arg_offset)) - k = args(1+arg_offset).nint_value (); - - if (!error_state && nargin > (2+arg_offset)) - if (args(2+arg_offset).is_real_scalar () || - args(2+arg_offset).is_complex_scalar ()) - { - sigma = args(2+arg_offset).complex_value (); - have_sigma = true; - } - else if (args(2+arg_offset).is_string ()) - { - typ = args(2+arg_offset).string_value (); - - // Use STL function to convert to upper case - transform (typ.begin (), typ.end (), typ.begin (), toupper); - - sigma = 0.; - } - else - error ("eigs: sigma must be a scalar or a string"); - - sigmar = std::real (sigma); - sigmai = std::imag (sigma); - - if (!error_state && nargin > (3+arg_offset)) - if (args(3+arg_offset).is_map ()) - { - Octave_map map = args(3+arg_offset).map_value (); - - // issym is ignored if A is not a function - if (map.contains ("issym")) - if (have_a_fun) - symmetric = map.contents ("issym")(0).double_value () != 0.; - - // isreal is ignored if A is not a function - if (map.contains ("isreal")) - if (have_a_fun) - a_is_complex = ! (map.contents ("isreal")(0).double_value () != 0.); - - if (map.contains ("tol")) - tol = map.contents ("tol")(0).double_value (); - - if (map.contains ("maxit")) - maxit = map.contents ("maxit")(0).nint_value (); - - if (map.contains ("p")) - p = map.contents ("p")(0).nint_value (); - - if (map.contains ("v0")) - { - if (a_is_complex || b_is_complex) - cresid = ComplexColumnVector - (map.contents ("v0")(0).complex_vector_value()); - else - resid = ColumnVector (map.contents ("v0")(0).vector_value()); - } - - if (map.contains ("disp")) - disp = map.contents ("disp")(0).nint_value (); - - if (map.contains ("cholB")) - cholB = map.contents ("cholB")(0).double_value () != 0.; - - if (map.contains ("permB")) - permB = ColumnVector (map.contents ("permB")(0).vector_value ()) - - 1.0; - } - else - error ("eigs: options argument must be a structure"); - - if (nargin > (4+arg_offset)) - error ("eigs: incorrect number of arguments"); - - if (have_b) - { - if (a_is_complex || b_is_complex) - { - if (a_is_sparse) - bscm = args(b_arg).sparse_complex_matrix_value (); - else - bcm = args(b_arg).complex_matrix_value (); - } - else - { - if (a_is_sparse) - bsmm = args(b_arg).sparse_matrix_value (); - else - bmm = args(b_arg).matrix_value (); - } - } - - // Mode 1 for SM mode seems unstable for some reason. - // Use Mode 3 instead, with sigma = 0. - if (!error_state && !have_sigma && typ == "SM") - have_sigma = true; - - if (!error_state) - { - octave_idx_type nconv; - if (a_is_complex || b_is_complex) - { - ComplexMatrix eig_vec; - ComplexColumnVector eig_val; - - - if (have_a_fun) - nconv = EigsComplexNonSymmetricFunc - (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, eig_val, - cresid, octave_stdout, tol, (nargout > 1), cholB, disp, maxit); - else if (have_sigma) - if (a_is_sparse) - nconv = EigsComplexNonSymmetricMatrixShift - (ascm, sigma, k, p, info, eig_vec, eig_val, bscm, permB, - cresid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - else - nconv = EigsComplexNonSymmetricMatrixShift - (acm, sigma, k, p, info, eig_vec, eig_val, bcm, permB, cresid, - octave_stdout, tol, (nargout > 1), cholB, disp, maxit); - else - if (a_is_sparse) - nconv = EigsComplexNonSymmetricMatrix - (ascm, typ, k, p, info, eig_vec, eig_val, bscm, permB, cresid, - octave_stdout, tol, (nargout > 1), cholB, disp, maxit); - else - nconv = EigsComplexNonSymmetricMatrix - (acm, typ, k, p, info, eig_vec, eig_val, bcm, permB, cresid, - octave_stdout, tol, (nargout > 1), cholB, disp, maxit); - - if (nargout < 2) - retval (0) = eig_val; - else - { - retval(2) = double (info); - retval(1) = ComplexDiagMatrix (eig_val); - retval(0) = eig_vec; - } - } - else if (sigmai != 0.) - { - // Promote real problem to a complex one. - ComplexMatrix eig_vec; - ComplexColumnVector eig_val; - - if (have_a_fun) - nconv = EigsComplexNonSymmetricFunc - (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec, eig_val, - cresid, octave_stdout, tol, (nargout > 1), cholB, disp, maxit); - else - if (a_is_sparse) - nconv = EigsComplexNonSymmetricMatrixShift - (SparseComplexMatrix (asmm), sigma, k, p, info, eig_vec, - eig_val, SparseComplexMatrix (bsmm), permB, cresid, - octave_stdout, tol, (nargout > 1), cholB, disp, maxit); - else - nconv = EigsComplexNonSymmetricMatrixShift - (ComplexMatrix (amm), sigma, k, p, info, eig_vec, - eig_val, ComplexMatrix (bmm), permB, cresid, - octave_stdout, tol, (nargout > 1), cholB, disp, maxit); - - if (nargout < 2) - retval (0) = eig_val; - else - { - retval(2) = double (info); - retval(1) = ComplexDiagMatrix (eig_val); - retval(0) = eig_vec; - } - } - else - { - if (symmetric) - { - Matrix eig_vec; - ColumnVector eig_val; - - if (have_a_fun) - nconv = EigsRealSymmetricFunc - (eigs_func, n, typ, sigmar, k, p, info, eig_vec, eig_val, - resid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - else if (have_sigma) - if (a_is_sparse) - nconv = EigsRealSymmetricMatrixShift - (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm, permB, - resid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - else - nconv = EigsRealSymmetricMatrixShift - (amm, sigmar, k, p, info, eig_vec, eig_val, bmm, permB, - resid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - else - if (a_is_sparse) - nconv = EigsRealSymmetricMatrix - (asmm, typ, k, p, info, eig_vec, eig_val, bsmm, permB, - resid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - else - nconv = EigsRealSymmetricMatrix - (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB, - resid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - - if (nargout < 2) - retval (0) = eig_val; - else - { - retval(2) = double (info); - retval(1) = DiagMatrix (eig_val); - retval(0) = eig_vec; - } - } - else - { - ComplexMatrix eig_vec; - ComplexColumnVector eig_val; - - if (have_a_fun) - nconv = EigsRealNonSymmetricFunc - (eigs_func, n, typ, sigmar, k, p, info, eig_vec, eig_val, - resid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - else if (have_sigma) - if (a_is_sparse) - nconv = EigsRealNonSymmetricMatrixShift - (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm, permB, - resid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - else - nconv = EigsRealNonSymmetricMatrixShift - (amm, sigmar, k, p, info, eig_vec, eig_val, bmm, permB, - resid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - else - if (a_is_sparse) - nconv = EigsRealNonSymmetricMatrix - (asmm, typ, k, p, info, eig_vec, eig_val, bsmm, permB, - resid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - else - nconv = EigsRealNonSymmetricMatrix - (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB, - resid, octave_stdout, tol, (nargout > 1), cholB, disp, - maxit); - - if (nargout < 2) - retval (0) = eig_val; - else - { - retval(2) = double (info); - retval(1) = ComplexDiagMatrix (eig_val); - retval(0) = eig_vec; - } - } - } - - if (nconv <= 0) - warning ("eigs: None of the %d requested eigenvalues converged", k); - else if (nconv < k) - warning ("eigs: Only %d of the %d requested eigenvalues converged", - nconv, k); - } - - if (! fcn_name.empty ()) - clear_function (fcn_name); - - unwind_protect::run_frame ("Feigs"); -#else - error ("eigs: not available in this version of Octave"); -#endif - - return retval; -} - -/* #### SPARSE MATRIX VERSIONS #### */ - -/* - -%% Real positive definite tests, n must be even -%!shared n, k, A, d0, d2 -%! n = 20; -%! k = 4; -%! A = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)]); -%! d0 = eig (A); -%! d2 = sort(d0); -%! [dum, idx] = sort( abs(d0)); -%! d0 = d0(idx); -%!test -%! d1 = eigs (A, k); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); -%!test -%! d1 = eigs (A,k+1); -%! assert (d1, d0(end:-1:(end-k)),1e-12); -%!test -%! d1 = eigs (A, k, 'lm'); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); -%!test -%! d1 = eigs (A, k, 'sm'); -%! assert (d1, d0(k:-1:1), 1e-12); -%!test -%! d1 = eigs (A, k, 'la'); -%! assert (d1, d2(end:-1:(end-k+1)), 1e-12); -%!test -%! d1 = eigs (A, k, 'sa'); -%! assert (d1, d2(1:k), 1e-12); -%!test -%! d1 = eigs (A, k, 'be'); -%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-12); -%!test -%! d1 = eigs (A, k+1, 'be'); -%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-12); -%!test -%! d1 = eigs (A, k, 4.1); -%! [dum,idx0] = sort (abs(d0 - 4.1)); -%! [dum,idx1] = sort (abs(d1 - 4.1)); -%! assert (d1(idx1), d0(idx0(1:k)), 1e-12); -%!test -%! d1 = eigs(A, speye(n), k, 'lm'); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!assert (eigs(A,k,4.1), eigs(A,speye(n),k,4.1), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, speye(n), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, speye(n), k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!assert (eigs(A,k,4.1), eigs(A,speye(n),k,4.1), 1e-12); -%!test -%! fn = @(x) A * x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 'lm', opts); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); -%!test -%! fn = @(x) A \ x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 'sm', opts); -%! assert (d1, d0(k:-1:1), 1e-12); -%!test -%! fn = @(x) (A - 4.1 * eye(n)) \ x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (d1, eigs(A,k,4.1), 1e-12); -%!test -%! [v1,d1] = eigs(A, k, 'lm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'la'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sa'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'be'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor - -*/ - -/* - -%% Real unsymmetric tests -%!shared n, k, A, d0 -%! n = 20; -%! k = 4; -%! A = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)]); -%! d0 = eig (A); -%! [dum, idx] = sort(abs(d0)); -%! d0 = d0(idx); -%!test -%! d1 = eigs (A, k); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A,k+1); -%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12); -%!test -%! d1 = eigs (A, k, 'lm'); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A, k, 'sm'); -%! assert (abs(d1), abs(d0(1:k)), 1e-12); -%!test -%! d1 = eigs (A, k, 'lr'); -%! [dum, idx] = sort (real(d0)); -%! d2 = d0(idx); -%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A, k, 'sr'); -%! [dum, idx] = sort (real(abs(d0))); -%! d2 = d0(idx); -%! assert (real(d1), real(d2(1:k)), 1e-12); -%!test -%! d1 = eigs (A, k, 'li'); -%! [dum, idx] = sort (imag(abs(d0))); -%! d2 = d0(idx); -%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12); -%!test -%! d1 = eigs (A, k, 'si'); -%! [dum, idx] = sort (imag(abs(d0))); -%! d2 = d0(idx); -%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12); -%!test -%! d1 = eigs (A, k, 4.1); -%! [dum,idx0] = sort (abs(d0 - 4.1)); -%! [dum,idx1] = sort (abs(d1 - 4.1)); -%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12); -%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12); -%!test -%! d1 = eigs(A, speye(n), k, 'lm'); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, speye(n), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, speye(n), k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,speye(n),k,4.1)), 1e-12); -%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,speye(n),k,4.1))), 1e-12); -%!test -%! fn = @(x) A * x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! fn = @(x) A \ x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 'sm', opts); -%! assert (abs(d1), d0(1:k), 1e-12); -%!test -%! fn = @(x) (A - 4.1 * eye(n)) \ x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!test -%! [v1,d1] = eigs(A, k, 'lm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'lr'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sr'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'li'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'si'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor - -*/ - -/* - -%% Complex hermitian tests -%!shared n, k, A, d0 -%! n = 20; -%! k = 4; -%! A = sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)]); -%! d0 = eig (A); -%! [dum, idx] = sort(abs(d0)); -%! d0 = d0(idx); -%!test -%! d1 = eigs (A, k); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A,k+1); -%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12); -%!test -%! d1 = eigs (A, k, 'lm'); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A, k, 'sm'); -%! assert (abs(d1), abs(d0(1:k)), 1e-12); -%!test -%! d1 = eigs (A, k, 'lr'); -%! [dum, idx] = sort (real(abs(d0))); -%! d2 = d0(idx); -%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A, k, 'sr'); -%! [dum, idx] = sort (real(abs(d0))); -%! d2 = d0(idx); -%! assert (real(d1), real(d2(1:k)), 1e-12); -%!test -%! d1 = eigs (A, k, 'li'); -%! [dum, idx] = sort (imag(abs(d0))); -%! d2 = d0(idx); -%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12); -%!test -%! d1 = eigs (A, k, 'si'); -%! [dum, idx] = sort (imag(abs(d0))); -%! d2 = d0(idx); -%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12); -%!test -%! d1 = eigs (A, k, 4.1); -%! [dum,idx0] = sort (abs(d0 - 4.1)); -%! [dum,idx1] = sort (abs(d1 - 4.1)); -%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12); -%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12); -%!test -%! d1 = eigs(A, speye(n), k, 'lm'); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, speye(n), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, speye(n)(q,q), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, speye(n), k, 4.1, opts); -%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12); -%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts); -%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12); -%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12); -%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,speye(n),k,4.1)), 1e-12); -%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,speye(n),k,4.1))), 1e-12); -%!test -%! fn = @(x) A * x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! fn = @(x) A \ x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, 'sm', opts); -%! assert (abs(d1), d0(1:k), 1e-12); -%!test -%! fn = @(x) (A - 4.1 * eye(n)) \ x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!test -%! [v1,d1] = eigs(A, k, 'lm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'lr'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sr'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'li'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'si'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*speye(n))*v1(:,i))),0.,1e-12) -%! endfor - -*/ - -/* #### FULL MATRIX VERSIONS #### */ - -/* - -%% Real positive definite tests, n must be even -%!shared n, k, A, d0, d2 -%! n = 20; -%! k = 4; -%! A = full(sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),4*ones(1,n),ones(1,n-2)])); -%! d0 = eig (A); -%! d2 = sort(d0); -%! [dum, idx] = sort( abs(d0)); -%! d0 = d0(idx); -%!test -%! d1 = eigs (A, k); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); -%!test -%! d1 = eigs (A,k+1); -%! assert (d1, d0(end:-1:(end-k)),1e-12); -%!test -%! d1 = eigs (A, k, 'lm'); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); -%!test -%! d1 = eigs (A, k, 'sm'); -%! assert (d1, d0(k:-1:1), 1e-12); -%!test -%! d1 = eigs (A, k, 'la'); -%! assert (d1, d2(end:-1:(end-k+1)), 1e-12); -%!test -%! d1 = eigs (A, k, 'sa'); -%! assert (d1, d2(1:k), 1e-12); -%!test -%! d1 = eigs (A, k, 'be'); -%! assert (d1, d2([1:floor(k/2), (end - ceil(k/2) + 1):end]), 1e-12); -%!test -%! d1 = eigs (A, k+1, 'be'); -%! assert (d1, d2([1:floor((k+1)/2), (end - ceil((k+1)/2) + 1):end]), 1e-12); -%!test -%! d1 = eigs (A, k, 4.1); -%! [dum,idx0] = sort (abs(d0 - 4.1)); -%! [dum,idx1] = sort (abs(d1 - 4.1)); -%! assert (d1(idx1), d0(idx0(1:k)), 1e-12); -%!test -%! d1 = eigs(A, eye(n), k, 'lm'); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!assert (eigs(A,k,4.1), eigs(A,eye(n),k,4.1), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, eye(n), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, eye(n), k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!assert (eigs(A,k,4.1), eigs(A,eye(n),k,4.1), 1e-12); -%!test -%! fn = @(x) A * x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 'lm', opts); -%! assert (d1, d0(end:-1:(end-k+1)), 1e-12); -%!test -%! fn = @(x) A \ x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 'sm', opts); -%! assert (d1, d0(k:-1:1), 1e-12); -%!test -%! fn = @(x) (A - 4.1 * eye(n)) \ x; -%! opts.issym = 1; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (d1, eigs(A,k,4.1), 1e-12); -%!test -%! [v1,d1] = eigs(A, k, 'lm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'la'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sa'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'be'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor - -*/ - -/* - -%% Real unsymmetric tests -%!shared n, k, A, d0 -%! n = 20; -%! k = 4; -%! A = full(sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[ones(1,n-2),1:n,-ones(1,n-2)])); -%! d0 = eig (A); -%! [dum, idx] = sort(abs(d0)); -%! d0 = d0(idx); -%!test -%! d1 = eigs (A, k); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A,k+1); -%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12); -%!test -%! d1 = eigs (A, k, 'lm'); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A, k, 'sm'); -%! assert (abs(d1), abs(d0(1:k)), 1e-12); -%!test -%! d1 = eigs (A, k, 'lr'); -%! [dum, idx] = sort (real(d0)); -%! d2 = d0(idx); -%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A, k, 'sr'); -%! [dum, idx] = sort (real(abs(d0))); -%! d2 = d0(idx); -%! assert (real(d1), real(d2(1:k)), 1e-12); -%!test -%! d1 = eigs (A, k, 'li'); -%! [dum, idx] = sort (imag(abs(d0))); -%! d2 = d0(idx); -%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12); -%!test -%! d1 = eigs (A, k, 'si'); -%! [dum, idx] = sort (imag(abs(d0))); -%! d2 = d0(idx); -%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12); -%!test -%! d1 = eigs (A, k, 4.1); -%! [dum,idx0] = sort (abs(d0 - 4.1)); -%! [dum,idx1] = sort (abs(d1 - 4.1)); -%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12); -%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12); -%!test -%! d1 = eigs(A, eye(n), k, 'lm'); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, eye(n), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, eye(n), k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,eye(n),k,4.1)), 1e-12); -%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,eye(n),k,4.1))), 1e-12); -%!test -%! fn = @(x) A * x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! fn = @(x) A \ x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 'sm', opts); -%! assert (abs(d1), d0(1:k), 1e-12); -%!test -%! fn = @(x) (A - 4.1 * eye(n)) \ x; -%! opts.issym = 0; opts.isreal = 1; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!test -%! [v1,d1] = eigs(A, k, 'lm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'lr'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sr'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'li'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'si'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor - -*/ - -/* - -%% Complex hermitian tests -%!shared n, k, A, d0 -%! n = 20; -%! k = 4; -%! A = full(sparse([3:n,1:n,1:(n-2)],[1:(n-2),1:n,3:n],[1i*ones(1,n-2),4*ones(1,n),-1i*ones(1,n-2)])); -%! d0 = eig (A); -%! [dum, idx] = sort(abs(d0)); -%! d0 = d0(idx); -%!test -%! d1 = eigs (A, k); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A,k+1); -%! assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12); -%!test -%! d1 = eigs (A, k, 'lm'); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A, k, 'sm'); -%! assert (abs(d1), abs(d0(1:k)), 1e-12); -%!test -%! d1 = eigs (A, k, 'lr'); -%! [dum, idx] = sort (real(abs(d0))); -%! d2 = d0(idx); -%! assert (real(d1), real(d2(end:-1:(end-k+1))), 1e-12); -%!test -%! d1 = eigs (A, k, 'sr'); -%! [dum, idx] = sort (real(abs(d0))); -%! d2 = d0(idx); -%! assert (real(d1), real(d2(1:k)), 1e-12); -%!test -%! d1 = eigs (A, k, 'li'); -%! [dum, idx] = sort (imag(abs(d0))); -%! d2 = d0(idx); -%! assert (sort(imag(d1)), sort(imag(d2(end:-1:(end-k+1)))), 1e-12); -%!test -%! d1 = eigs (A, k, 'si'); -%! [dum, idx] = sort (imag(abs(d0))); -%! d2 = d0(idx); -%! assert (sort(imag(d1)), sort(imag(d2(1:k))), 1e-12); -%!test -%! d1 = eigs (A, k, 4.1); -%! [dum,idx0] = sort (abs(d0 - 4.1)); -%! [dum,idx1] = sort (abs(d1 - 4.1)); -%! assert (abs(d1(idx1)), abs(d0(idx0(1:k))), 1e-12); -%! assert (sort(imag(d1(idx1))), sort(imag(d0(idx0(1:k)))), 1e-12); -%!test -%! d1 = eigs(A, eye(n), k, 'lm'); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, eye(n), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, eye(n)(q,q), k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! opts.cholB=true; -%! d1 = eigs(A, eye(n), k, 4.1, opts); -%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12); -%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12); -%!test -%! opts.cholB=true; -%! q = [2:n,1]; -%! opts.permB=q; -%! d1 = eigs(A, eye(n)(q,q), k, 4.1, opts); -%! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12); -%! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12); -%!assert (abs(eigs(A,k,4.1)), abs(eigs(A,eye(n),k,4.1)), 1e-12); -%!assert (sort(imag(eigs(A,k,4.1))), sort(imag(eigs(A,eye(n),k,4.1))), 1e-12); -%!test -%! fn = @(x) A * x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, 'lm', opts); -%! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12); -%!test -%! fn = @(x) A \ x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, 'sm', opts); -%! assert (abs(d1), d0(1:k), 1e-12); -%!test -%! fn = @(x) (A - 4.1 * eye(n)) \ x; -%! opts.issym = 0; opts.isreal = 0; -%! d1 = eigs (fn, n, k, 4.1, opts); -%! assert (abs(d1), eigs(A,k,4.1), 1e-12); -%!test -%! [v1,d1] = eigs(A, k, 'lm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sm'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'lr'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'sr'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'li'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor -%!test -%! [v1,d1] = eigs(A, k, 'si'); -%! d1 = diag(d1); -%! for i=1:k -%! assert(max(abs((A - d1(i)*eye(n))*v1(:,i))),0.,1e-12) -%! endfor - -*/ - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/