 One web page per Section   All as one web page Math in Latex   f2cad.htm _printable.htm Math in MathML   f2cad.xml _priintable.xml

a: Purpose

a.b: Other Libraries and C++ AD Types
Provide the infrastructure to convert any Fortran 77 code to any C++ AD type.

b: get_started
The getting started examples in 4: get_started show how to use the files in the library with each of the AD types.

c: Example
The examples and tests in 6: library demonstrate the use of each of the Fortran routines in the library with each of the AD types. On the other hand, you should use the interface best for your particular AD type when using the library (see 4: get_started .)

d: Contents

f2cad-20100428: Convert Fortran 77 to C++ AD Types: : f2cad
Getting Started Using the f2cad Libraries: 4: get_started
C++ Types Corresponding to Fortran Prototypes: 5: prototype
d9knus: 6.1: d9knus
d9knus.f: 6.1.1: d9knus.f
d9lgmc: 6.2: d9lgmc
d9lgmc.f: 6.2.1: d9lgmc.f
daxpy: 6.3: daxpy
daxpy.f: 6.3.1: daxpy.f
dbesks: 6.4: dbesks
dbesks.f: 6.4.1: dbesks.f
dbskes: 6.5: dbskes
dbskes.f: 6.5.1: dbskes.f
dcsevl: 6.6: dcsevl
dcsevl.f: 6.6.1: dcsevl.f
ddot: 6.7: ddot
ddot.f: 6.7.1: ddot.f
dgamlm: 6.8: dgamlm
dgamlm.f: 6.8.1: dgamlm.f
dgamma: 6.9: dgamma
dgamma.f: 6.9.1: dgamma.f
dgefa: 6.10: dgefa
dgefa.f: 6.10.1: dgefa.f
dgesl: 6.11: dgesl
dgesl.f: 6.11.1: dgesl.f
dint: 6.12: dint
dint.f: 6.12.1: dint.f
dscal: 6.13: dscal
dscal.f: 6.13.1: dscal.f
idamax: 6.14: idamax
idamax.f: 6.14.1: idamax.f
initds: 6.15: initds
initds.f: 6.15.1: initds.f
xerror: 6.16: xerror
xerror.f: 6.16.1: xerror.f
Run Test of All the Library Routines: 7: run.cpp
Utilities used to Test Library: 8: utility
Implementation of Cmlib i1mach Routine in Cpp: 8.1: i1mach
Implementation of Cmlib d1mach Routine in Cpp: 8.2: d1mach
Determine if Two Double Values Are Nearly Equal: 8.4: near_equal
Implementation of Cmlib xerror Routine in Cpp: 8.5: xerror2cpp
Alphabetic Listing of Cross Reference Tags: 11: _reference
Keyword Index: 12: _index
External Internet References: 13: _external


This software is licensed to you under the terms of the GNU General Public License. A brief description of the license follows:  Below is a copy of the GNU General Public License that is provided along with this software. It is a verbatim copy of the file        http://www.gnu.org/licenses/gpl.txt

2.b: Complete Software License  GNU GENERAL PUBLIC LICENSE Version 2, June 1991 Copyright (C) 1989, 1991 Free Software Foundation, Inc. 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Everyone is permitted to copy and distribute verbatim copies of this license document, but changing it is not allowed. Preamble The licenses for most software are designed to take away your freedom to share and change it. By contrast, the GNU General Public License is intended to guarantee your freedom to share and change free software--to make sure the software is free for all its users. This General Public License applies to most of the Free Software Foundation's software and to any other program whose authors commit to using it. (Some other Free Software Foundation software is covered by the GNU Library General Public License instead.) You can apply it to your programs, too. When we speak of free software, we are referring to freedom, not price. Our General Public Licenses are designed to make sure that you have the freedom to distribute copies of free software (and charge for this service if you wish), that you receive source code or can get it if you want it, that you can change the software or use pieces of it in new free programs; and that you know you can do these things. To protect your rights, we need to make restrictions that forbid anyone to deny you these rights or to ask you to surrender the rights. These restrictions translate to certain responsibilities for you if you distribute copies of the software, or if you modify it. For example, if you distribute copies of such a program, whether gratis or for a fee, you must give the recipients all the rights that you have. You must make sure that they, too, receive or can get the source code. And you must show them these terms so they know their rights. We protect your rights with two steps: (1) copyright the software, and (2) offer you this license which gives you legal permission to copy, distribute and/or modify the software. Also, for each author's protection and ours, we want to make certain that everyone understands that there is no warranty for this free software. If the software is modified by someone else and passed on, we want its recipients to know that what they have is not the original, so that any problems introduced by others will not reflect on the original authors' reputations. Finally, any free program is threatened constantly by software patents. We wish to avoid the danger that redistributors of a free program will individually obtain patent licenses, in effect making the program proprietary. To prevent this, we have made it clear that any patent must be licensed for everyone's free use or not licensed at all. The precise terms and conditions for copying, distribution and modification follow. GNU GENERAL PUBLIC LICENSE TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION 0. This License applies to any program or other work which contains a notice placed by the copyright holder saying it may be distributed under the terms of this General Public License. The "Program", below, refers to any such program or work, and a "work based on the Program" means either the Program or any derivative work under copyright law: that is to say, a work containing the Program or a portion of it, either verbatim or with modifications and/or translated into another language. (Hereinafter, translation is included without limitation in the term "modification".) Each licensee is addressed as "you". Activities other than copying, distribution and modification are not covered by this License; they are outside its scope. The act of running the Program is not restricted, and the output from the Program is covered only if its contents constitute a work based on the Program (independent of having been made by running the Program). Whether that is true depends on what the Program does. 1. You may copy and distribute verbatim copies of the Program's source code as you receive it, in any medium, provided that you conspicuously and appropriately publish on each copy an appropriate copyright notice and disclaimer of warranty; keep intact all the notices that refer to this License and to the absence of any warranty; and give any other recipients of the Program a copy of this License along with the Program. You may charge a fee for the physical act of transferring a copy, and you may at your option offer warranty protection in exchange for a fee. 2. You may modify your copy or copies of the Program or any portion of it, thus forming a work based on the Program, and copy and distribute such modifications or work under the terms of Section 1 above, provided that you also meet all of these conditions: a) You must cause the modified files to carry prominent notices stating that you changed the files and the date of any change. b) You must cause any work that you distribute or publish, that in whole or in part contains or is derived from the Program or any part thereof, to be licensed as a whole at no charge to all third parties under the terms of this License. c) If the modified program normally reads commands interactively when run, you must cause it, when started running for such interactive use in the most ordinary way, to print or display an announcement including an appropriate copyright notice and a notice that there is no warranty (or else, saying that you provide a warranty) and that users may redistribute the program under these conditions, and telling the user how to view a copy of this License. (Exception: if the Program itself is interactive but does not normally print such an announcement, your work based on the Program is not required to print an announcement.) These requirements apply to the modified work as a whole. If identifiable sections of that work are not derived from the Program, and can be reasonably considered independent and separate works in themselves, then this License, and its terms, do not apply to those sections when you distribute them as separate works. But when you distribute the same sections as part of a whole which is a work based on the Program, the distribution of the whole must be on the terms of this License, whose permissions for other licensees extend to the entire whole, and thus to each and every part regardless of who wrote it. Thus, it is not the intent of this section to claim rights or contest your rights to work written entirely by you; rather, the intent is to exercise the right to control the distribution of derivative or collective works based on the Program. In addition, mere aggregation of another work not based on the Program with the Program (or with a work based on the Program) on a volume of a storage or distribution medium does not bring the other work under the scope of this License. 3. You may copy and distribute the Program (or a work based on it, under Section 2) in object code or executable form under the terms of Sections 1 and 2 above provided that you also do one of the following: a) Accompany it with the complete corresponding machine-readable source code, which must be distributed under the terms of Sections 1 and 2 above on a medium customarily used for software interchange; or, b) Accompany it with a written offer, valid for at least three years, to give any third party, for a charge no more than your cost of physically performing source distribution, a complete machine-readable copy of the corresponding source code, to be distributed under the terms of Sections 1 and 2 above on a medium customarily used for software interchange; or, c) Accompany it with the information you received as to the offer to distribute corresponding source code. (This alternative is allowed only for noncommercial distribution and only if you received the program in object code or executable form with such an offer, in accord with Subsection b above.) The source code for a work means the preferred form of the work for making modifications to it. For an executable work, complete source code means all the source code for all modules it contains, plus any associated interface definition files, plus the scripts used to control compilation and installation of the executable. However, as a special exception, the source code distributed need not include anything that is normally distributed (in either source or binary form) with the major components (compiler, kernel, and so on) of the operating system on which the executable runs, unless that component itself accompanies the executable. If distribution of executable or object code is made by offering access to copy from a designated place, then offering equivalent access to copy the source code from the same place counts as distribution of the source code, even though third parties are not compelled to copy the source along with the object code. 4. You may not copy, modify, sublicense, or distribute the Program except as expressly provided under this License. Any attempt otherwise to copy, modify, sublicense or distribute the Program is void, and will automatically terminate your rights under this License. However, parties who have received copies, or rights, from you under this License will not have their licenses terminated so long as such parties remain in full compliance. 5. You are not required to accept this License, since you have not signed it. However, nothing else grants you permission to modify or distribute the Program or its derivative works. These actions are prohibited by law if you do not accept this License. Therefore, by modifying or distributing the Program (or any work based on the Program), you indicate your acceptance of this License to do so, and all its terms and conditions for copying, distributing or modifying the Program or works based on it. 6. Each time you redistribute the Program (or any work based on the Program), the recipient automatically receives a license from the original licensor to copy, distribute or modify the Program subject to these terms and conditions. You may not impose any further restrictions on the recipients' exercise of the rights granted herein. You are not responsible for enforcing compliance by third parties to this License. 7. If, as a consequence of a court judgment or allegation of patent infringement or for any other reason (not limited to patent issues), conditions are imposed on you (whether by court order, agreement or otherwise) that contradict the conditions of this License, they do not excuse you from the conditions of this License. If you cannot distribute so as to satisfy simultaneously your obligations under this License and any other pertinent obligations, then as a consequence you may not distribute the Program at all. For example, if a patent license would not permit royalty-free redistribution of the Program by all those who receive copies directly or indirectly through you, then the only way you could satisfy both it and this License would be to refrain entirely from distribution of the Program. If any portion of this section is held invalid or unenforceable under any particular circumstance, the balance of the section is intended to apply and the section as a whole is intended to apply in other circumstances. It is not the purpose of this section to induce you to infringe any patents or other property right claims or to contest validity of any such claims; this section has the sole purpose of protecting the integrity of the free software distribution system, which is implemented by public license practices. Many people have made generous contributions to the wide range of software distributed through that system in reliance on consistent application of that system; it is up to the author/donor to decide if he or she is willing to distribute software through any other system and a licensee cannot impose that choice. This section is intended to make thoroughly clear what is believed to be a consequence of the rest of this License. 8. If the distribution and/or use of the Program is restricted in certain countries either by patents or by copyrighted interfaces, the original copyright holder who places the Program under this License may add an explicit geographical distribution limitation excluding those countries, so that distribution is permitted only in or among countries not thus excluded. In such case, this License incorporates the limitation as if written in the body of this License. 9. The Free Software Foundation may publish revised and/or new versions of the General Public License from time to time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or concerns. Each version is given a distinguishing version number. If the Program specifies a version number of this License which applies to it and "any later version", you have the option of following the terms and conditions either of that version or of any later version published by the Free Software Foundation. If the Program does not specify a version number of this License, you may choose any version ever published by the Free Software Foundation. 10. If you wish to incorporate parts of the Program into other free programs whose distribution conditions are different, write to the author to ask for permission. For software which is copyrighted by the Free Software Foundation, write to the Free Software Foundation; we sometimes make exceptions for this. Our decision will be guided by the two goals of preserving the free status of all derivatives of our free software and of promoting the sharing and reuse of software generally. NO WARRANTY 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. END OF TERMS AND CONDITIONS How to Apply These Terms to Your New Programs If you develop a new program, and you want it to be of the greatest possible use to the public, the best way to achieve this is to make it free software which everyone can redistribute and change under these terms. To do so, attach the following notices to the program. It is safest to attach them to the start of each source file to most effectively convey the exclusion of warranty; and each file should have at least the "copyright" line and a pointer to where the full notice is found. <one line to give the program's name and a brief idea of what it does.> Copyright (C) <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, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Also add information on how to contact you by electronic and paper mail. If the program is interactive, make it output a short notice like this when it starts in an interactive mode: Gnomovision version 69, Copyright (C) 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. 

3.a: Subversion
If you have subversion (http://subversion.tigris.org/) , you can use the following command to get the specified release of       svn checkout $f2cad_svn/releases/0.yyyymmdd.r  where yyyy is the year, mm is the month, dd is the day, and r is the release number. 3.b: Tarballs 3.b.a: Download You can download a tarball release of the form f2cAD-0.yyyymmdd.r (http://www.coin-or.org/download/source/CoinBazaar/) , 3.b.b: Extraction Use the command  tar -xvzf f2cAD-0.yyyymmdd.r.tgz  to extract the download file into the distribution directory  f2cAD-0.yyyymmdd.r  3.c: User Documentation The user documentation for releases is included in the doc sub-directory of the distribution directory. The file doc/doc.xml is the top of the documentation tree (note table at top listing the different versions of the documentation). 3.d: Built Sources The next step is to use f2c (http://netlib.sandia.gov/f2c/index.html) to automatically build the C++ source code from the Fortran source code. In the distribution directory, execute the command:  ./build_sources.sh  3.e: Configure Enter the directory created by the extraction and execute the command:  ./configure \ [ prefix=prefix ] \ F2C_DIR=f2c_dir \ [ ADOLC_DIR=adolc_dir ] \ [ CPPAD_DIR=cppad_dir ] \ [ FADBAD_DIR=fadbad_dir ] \ [ CXX_FLAGS=cxx_flags ]  where values in italic are replaced values that you choose, and entries of the form [ ... ] are optional. 3.f: Prefix The default value for prefix is $HOME; i.e., by default the f2cad include files will install in       $HOME/include/f2cad  and the f2cad libraries are $HOME/lib/libf2c_adolc.a      $HOME/lib/libf2c_cppad.a$HOME/lib/libf2c_fadbad.a  If you want to install elsewhere, you will have to use this option. For example, if you have permission to write into the directories /usr/local/include and /usr/local/lib,  ./configure --prefix=/usr/local  will set up for installing the f2cad include files in the directory       /usr/local/include/f2cad  and the f2cad libraries will be       /usr/local/lib/libf2c_adolc.a      /usr/local/lib/libf2c_cppad.a      /usr/local/lib/libf2c_fadbad.a  This is the standard location for such files on many systems. If prefix is not a standard location for your system, you will need the compiler switches       -I prefix/include -L prefix/lib  on the command line when you compile and link your programs with the f2cad libraries.

3.g: F2C_DIR
The directory f2c_dir must be the location where f2c (http://netlib.sandia.gov/f2c/index.html) is installed. To be specific, the files       f2c_dir/include/f2c.h      f2c_dir/lib/libf2c.a      f2c_dir/bin/f2c  must be the include file, library file, and executable for f2c.

If you have adolc (https://projects.coin-or.org/ADOL-C) installed on your system, you can specify a value for adolc_dir in the 3.e: configure command line. The value of adolc_dir must be such that       adolc_dir/include/adolc/adouble.h  is a valid way to reference the include file adouble.h. In this case, the lib/adolc/libf2c_adolc.a library and the corresponding example and test program example/adolc/run will be built.

If you have cppad (http://www.coin-or.org/CppAD/) installed on your system, you can specify a value for cppad_dir in the 3.e: configure command line. The value of cppad_dir must be such that       cppad_dir/include/cppad/cppad.hpp  is a valid way to reference the include file cppad.hpp. In this case, the lib/cppad/libf2c_cppad.a library and the corresponding example and test program example/cppad/run will be built.

If you have fadbad (http://www.fadbad.com/) installed on your system, you can specify a value for fadbad_dir . It must be such that       fadbad_dir/FADBAD++/fadiff.h  is a valid reference to fadiff.h. In this case, the lib/fadbad/libf2c_fadbad.a library and the corresponding example and test program example/fadbad/run will be built.

3.k: C++ Errors And Warnings
If the command line argument cxx_flags is present, it specifies compiler error and warning flags. For example,       CXX_FLAGS="-Wall -ansi"  would specify that error and warning flags -Wall and -ansi should be included in all the C++ compile commands. The error and warning flags chosen must be valid options for the C++ compiler. The default value for cxx_flags is the empty string.

3.l: Building The Libraries
In the f2cad distribution directory, execute the command  make  This will build the libraries lib/adolc/libf2c_adolc.a, lib/cppad/libf2c_cppad.a, and lib/fadbad/libf2c_fadbad.a.

3.m: Testing
In the f2cad distribution directory, execute the command  make test  This will compile and run the tests. The results of the tests are stored in the file test.log.

3.n: Final Installation
Once you are satisfied that the tests are giving correct results, you can install the libraries and include files into easy to use directories by executing the command  make install  This will install f2cad in the location specified by 3.f: prefix .
Input File: omh/install.omh
4: Getting Started Using the f2cad Libraries

4.a: Description
The examples demonstrate how to use each of f2cad libraries libf2c_adolc.a, libf2c_cppad.a, and libf2c_fadbad.a. Adding to the libraries is more complicated and is described in the 9: add2lib.sh section.

4.b: Contents

Input File: omh/get_started.omh

4.1.a: Prototype
int f2cad::dscal_(integer *n, doublereal *da, doublereal *dx, integer *incx);

4.1.b: Description
This is a simple example using Adolc (http://www.math.tu-dresden.de/~adol-c/) and the libf2cadAdolc.a library routine 6.13: dscal to compute the function   $f(a) = a \left( \begin{array}{c} 1 \\ 2 \end{array} \right)$  We check that the derivative of this function, calculated using the Adolc, satisfies   $f^{(1)} (a) = \left( \begin{array}{c} 1 \\ 2 \end{array} \right)$   // standard input output library # include <iostream> // f2cad set up for adouble # define F2CAD_USE_ADOLC // prototype for dscal # include <f2cad/dscal.hpp> int main(void) { // turn on recording of adouble operations int tag = 1; int keep = 1; trace_on(tag, keep); // independent variables int n = 1; // number of independent variables adouble a[1]; // vector of independent variables a[0] <<= 5.; // first independent variables // create data structure expected by dscal integer m = 2; adouble f[2]; f[0] = 1.; f[1] = 2.; integer incf = 1; // set f = a * f f2cad::dscal_(&m, a, f, &incf ); // dependent variables double fi; f[0] >>= fi; // first dependent variable value f[1] >>= fi; // second dependent variable value // turn off recording of adouble operations trace_off(); // set differential for a int d = 1; // order of the differential double *daMemory = new double [2 * n]; double **da = new double * [n]; int i, j; for(j = 0; j < n; j++) da[j] = daMemory + 2 * j; da[0][0] = .5; // value of a da[0][1] = 1.; // differential for a // compute differential for f double *dfMemory = new double [2 * m]; double **df = new double * [m]; for(i = 0; i < m; i++) df[i] = dfMemory + 2 * i; forward(tag, int(m), n, d, keep, da, df); // print results std::cout << "f(a) = a * (1, 2)^T" << std::endl; std::cout << "f'(a) = (" << df[0][1] << ", " << df[1][1] << ")^T" << std::endl; bool ok = df[0][1] == 1. && df[1][1] == 2.; // free allocated memory delete [] daMemory; delete [] dfMemory; delete [] da; delete [] df; if( ok ) return 0; // no error return 1; // error } 

4.2.a: Prototype
int f2cad::dscal_(integer *n, doublereal *da, doublereal *dx, integer *incx);

4.2.b: Description
This is a simple example using CppAD (http://www.seanet.com/~bradbell/CppAD/) and the libf2cadCppAD.a library routine 6.13: dscal to compute the function   $f(a) = a \left( \begin{array}{c} 1 \\ 2 \end{array} \right)$  We check that the derivative of this function, calculated using the CppAD, satisfies   $f^{(1)} (a) = \left( \begin{array}{c} 1 \\ 2 \end{array} \right)$   // standard input output library # include <iostream> // f2cad set up for CppAD::AD<double> # define F2CAD_USE_CPPAD // prototype for dscal # include <f2cad/dscal.hpp> int main(void) { using namespace CppAD; // independent variables size_t n = 1; // number of independent variables CppADvector< AD<double> > a(n); // vector of independent variables a[0] = 5.; // value of independent variables Independent(a); // declare independent variables // create data structure expected by dscal AD<double> A[1]; A[0] = a[0]; integer M = 2; AD<double> F[2]; F[0] = 1.; F[1] = 2.; integer Incf = 1; // set f = a * f f2cad::dscal_(&M, A, F, &Incf ); // dependent variables size_t m = 2; // number of dependent variables CppADvector< AD<double> > f(2); // vector of dependent variables f[0] = F[0]; // value of dependent variables f[1] = F[1]; ADFun<double> Fun(a, f); // declare dependent variables // set differential for a CppADvector<double> da(n); CppADvector<double> df(m); da[0] = 1.; // compute differential for f df = Fun.Forward(1, da); // print results std::cout << "f(a) = a * (1, 2)^T" << std::endl; std::cout << "f'(a) = (" << df[0] << ", " << df[1] << ")^T" << std::endl; if( df[0] != 1. || df[1] != 2. ) return 1; // error return 0; // no error } 

4.3.a: Syntax
int f2cad::dscal_(integer *n, doublereal *da, doublereal *dx, integer *incx);

4.3.b: Description
This is a simple example using Fadbad (http://www.imm.dtu.dk/fadbad.html/) and the libf2cadFadbad.a library routine 6.13: dscal to compute the function   $f(a) = a \left( \begin{array}{c} 1 \\ 2 \end{array} \right)$  We check that the derivative of this function, calculated using the Fadbad, satisfies   $f^{(1)} (a) = \left( \begin{array}{c} 1 \\ 2 \end{array} \right)$   // standard input output library # include <iostream> // f2cad set up for F<double> # define F2CAD_USE_FADBAD // prototype for dscal # include <f2cad/dscal.hpp> int main(void) { // independent variables doublereal a[1]; // vector of independent variables a[0] = 5.; // value of independent variables a[0].diff(0, 1); // declare independent variables // create data structure expected by dscal integer m = 2; doublereal f[2]; f[0] = 1.; f[1] = 2.; integer incf = 1; // set f = a * f f2cad::dscal_(&m, a, f, &incf ); // print results std::cout << "f(a) = a * (1, 2)^T" << std::endl; std::cout << "f'(a) = (" << f[0].d(0) << ", " << f[1].d(0) << ")^T" << std::endl; // return value not significant if( f[0].d(0) != 1. || f[1].d(0) != 2. ) return 1; // error return 0; // no error } 
5: C++ Types Corresponding to Fortran Prototypes

5.a: Description
The following code (extracted from <f2cad/f2cad.hpp>) shows how the prototypes for Fortran functions link to C++ and the AD types:  typedef enum { test_none, test_pass, test_fail } test_result; typedef long int integer; # ifdef F2CAD_USE_ADOLC # include <adolc/adouble.h> # include <adolc/interfaces.h> # include <adolc/taping.h> typedef adouble doublereal; # endif # ifdef F2CAD_USE_CPPAD # include <cppad/cppad.hpp> typedef CppAD::AD<double> doublereal; # endif # ifdef F2CAD_USE_FADBAD # include <FADBAD++/fadiff.h> typedef fadbad::F<double> doublereal; # endif 
Input File: omh/prototype.omh

6.a: Description
This section contains examples and tests for each of the files in the library. These tests run with all three AD packages Adolc, CppAD, and Fadbad. For this reason they use each of the AD packages indirectly though the 8.3: f2cad_link interface In order to do so, these examples and tests include the file <f2cad/f2cad.h> and are not efficient. When you use the library in actual applications, you should not include <f2cad/f2cad.h> and instead define integer and doublereal directly (as done in the 4: get_started examples).

6.b: Contents

Input File: omh/library.omh
6.1: d9knus

6.1.a: Prototype
int f2cad::d9knus_(doublereal *xnu, doublereal *x, doublereal *bknu, doublereal *bknu1, integer *iswtch);

6.1.b: Fortran Source
6.1.1: d9knus.f

6.1.c: Notation
We use   K_{\nu} (x) to denote the modified Bessel function of the second kind; i.e., it is the solution of the ODE   $x^2 y^{(2)} (x) + x y^{(1)} (x) - ( x^2 + \nu^2 ) = 0$  that goes to zero as   x \rightarrow \infty .

6.1.d: Description
The routine d9knus computes the values   \exp(x) K_\nu (x) and   \exp(x) K_{\nu+1} (x) for   x \in [0,1] and   \nu \in [0,1] .

6.1.e: Special Case
By Formula 10.2.17 of Abramowitz and Stegun   $K_{1/2} ( x ) = \sqrt{ \pi / ( 2 x ) } \exp( - x )$   # include <cmath> # include <f2cad/d9knus.hpp> test_result d9knus(void) { bool ok = true; // independent variable integer n = 1; // number of independent variables doublereal x[1]; // vector of independent variables x[0] = .75; // value of independent variable f2cad::Independent(n, x); // declare x as independent variable vector // call the routine doublereal xnu = .5; doublereal bknu, bknu1; integer iswtch; f2cad::d9knus_(&xnu, x+0, &bknu, &bknu1, &iswtch); // dependent variables integer m = 1; // number of dependent variables doublereal f[1]; // vector of dependent variables f[0] = bknu; // value f[0] = K_nu (x[0]) * exp( x[0] ) f2cad::Dependent(m, f); // declare f as dependent variable vector // check function value double x0 = f2cad::Value( x[0] ); double f0 = f2cad::Value( f[0] ); double pi = 4. * std::atan(1.); double check = 1. / sqrt( 2. * x0 / pi ); ok &= f2cad::near_equal(f0, check, 1e-10, 1e-10); // Evaluate the derivative of f[0] w.r.t x[0] double p = f2cad::Partial<doublereal>(0, 0); check = - check * check * check / pi; ok &= f2cad::near_equal(p, check, 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/d9knus.cpp
6.1.1: d9knus.f
 subroutine d9knus(xnu,x,bknu,bknu1,iswtch) c***begin prologue d9knus c***date written 770601 (yymmdd) c***revision date 820801 (yymmdd) c***category no. c10b3 c***keywords bessel function,double precision,fnlib,special function c***author fullerton, w., (lanl) c***purpose computes bessel functions exp(x)*k-sub-xnu(x) and exp(x)* c k-sub-xnu+1(x) for 0. .le. xnu .lt. 1. c***description c c compute bessel functions exp(x) * k-sub-xnu (x) and c exp(x) * k-sub-xnu+1 (x) for 0.0 .le. xnu .lt. 1.0 . c c series for c0k on the interval 0. to 2.50000e-01 c with weighted error 2.16e-32 c log weighted error 31.67 c significant figures required 30.86 c decimal places required 32.40 c c series for znu1 on the interval -7.00000e-01 to 0. c with weighted error 2.45e-33 c log weighted error 32.61 c significant figures required 31.85 c decimal places required 33.26 c***references (none) c***routines called d1mach,dcsevl,dgamma,initds,xerror c***end prologue d9knus double precision xnu, x, bknu, bknu1, alpha(32), beta(32), a(32), 1 c0kcs(29), znu1cs(20), alnz, aln2, a0, bknud, bknu0, 2 b0, c0, euler, expx, p1, p2, p3, qq, result, sqpi2, sqrtx, v, 3 vlnz, xi, xmu, xnusml, xsml, x2n, x2tov, z, ztov, alnsml, 4 alnbig real alneps double precision d1mach, dcsevl, dgamma data c0k cs( 1) / +.6018305724 2626108387 5774451803 29 d-1 / data c0k cs( 2) / -.1536487143 3017286092 9597559431 24 d+0 / data c0k cs( 3) / -.1175117600 8210492040 0682292262 13 d-1 / data c0k cs( 4) / -.8524878889 1979509827 0484015509 87 d-3 / data c0k cs( 5) / -.6132983876 7496791874 0981769221 11 d-4 / data c0k cs( 6) / -.4405228124 5510444562 6798895485 05 d-5 / data c0k cs( 7) / -.3163124672 8384488192 9154458921 99 d-6 / data c0k cs( 8) / -.2271071938 2899588330 6737717933 96 d-7 / data c0k cs( 9) / -.1630564460 8077609552 2746205153 60 d-8 / data c0k cs( 10) / -.1170693929 9414776568 7560440431 30 d-9 / data c0k cs( 11) / -.8405206378 6464437174 5465934137 92 d-11 / data c0k cs( 12) / -.6034667011 8979991487 0960507371 98 d-12 / data c0k cs( 13) / -.4332696033 5681371952 0459973669 03 d-13 / data c0k cs( 14) / -.3110735803 0203546214 6346977722 37 d-14 / data c0k cs( 15) / -.2233407822 6736982254 4861334098 40 d-15 / data c0k cs( 16) / -.1603514671 6864226300 6357915286 10 d-16 / data c0k cs( 17) / -.1151271736 3666556196 0356977053 05 d-17 / data c0k cs( 18) / -.8265759174 6836959105 1694790892 58 d-19 / data c0k cs( 19) / -.5934548080 6383948172 3334366959 84 d-20 / data c0k cs( 20) / -.4260813819 6467143926 4996130239 76 d-21 / data c0k cs( 21) / -.3059126686 4812876299 2636983705 42 d-22 / data c0k cs( 22) / -.2196354142 6734575224 9755018155 16 d-23 / data c0k cs( 23) / -.1576911326 1495836071 1057506847 60 d-24 / data c0k cs( 24) / -.1132171393 5950320948 7577310480 56 d-25 / data c0k cs( 25) / -.8128624883 4598404082 7923497144 33 d-27 / data c0k cs( 26) / -.5836090089 3453226552 8293493159 49 d-28 / data c0k cs( 27) / -.4190124162 3610922519 4523377809 05 d-29 / data c0k cs( 28) / -.3008373796 0206435069 5305042128 62 d-30 / data c0k cs( 29) / -.2159915206 7808647728 3421680898 32 d-31 / data znu1cs( 1) / +.2033067569 9419172967 4444001216 911 d+0 / data znu1cs( 2) / +.1400779334 1321977106 2943670790 563 d+0 / data znu1cs( 3) / +.7916796961 0016135284 0972241972 320 d-2 / data znu1cs( 4) / +.3398011825 3210404535 2930092205 750 d-3 / data znu1cs( 5) / +.1174197568 8989336666 4507228352 690 d-4 / data znu1cs( 6) / +.3393575706 1226168033 3825865475 121 d-6 / data znu1cs( 7) / +.8425941769 7621991019 4629891264 803 d-8 / data znu1cs( 8) / +.1833366770 2485008918 4748150900 090 d-9 / data znu1cs( 9) / +.3549698447 0441631086 3007064469 557 d-11 / data znu1cs( 10) / +.6190324964 6988733220 5244342078 407 d-13 / data znu1cs( 11) / +.9819645356 8043942496 0346115456 527 d-15 / data znu1cs( 12) / +.1428513143 9649047421 1473563005 985 d-16 / data znu1cs( 13) / +.1918949218 8782529896 6162467488 436 d-18 / data znu1cs( 14) / +.2394309797 3949891416 2313140597 128 d-20 / data znu1cs( 15) / +.2788902468 1534735483 5870465474 995 d-22 / data znu1cs( 16) / +.3046066506 3303344258 2845214092 865 d-24 / data znu1cs( 17) / +.3131732370 4219181577 1564260932 089 d-26 / data znu1cs( 18) / +.3041330989 8785495164 5174908005 034 d-28 / data znu1cs( 19) / +.2798403846 3683308434 3185097659 733 d-30 / data znu1cs( 20) / +.2446371862 7449759648 5238794922 666 d-32 / data euler / 0.5772156649 0153286060 6512090082 40d0 / data sqpi2 / +1.253314137 3155002512 0788264240 55 d0 / data aln2 / 0.6931471805 5994530941 7232121458 18d0 / data ntc0k, ntznu1 / 2*0 / data xnusml, xsml, alnsml, alnbig / 4*0.d0 / data alneps / 0.0 / c***first executable statement d9knus if (ntc0k.ne.0) go to 10 eta = 0.1d0*d1mach(3) ntc0k = initds (c0kcs, 29, eta) ntznu1 = initds (znu1cs, 20, eta) c xnusml = dsqrt (d1mach(3)/8.d0) xsml = 0.1d0*d1mach(3) alnsml = dlog (d1mach(1)) alnbig = dlog (d1mach(2)) alneps = dlog (0.1d0*d1mach(3)) c 10 if (xnu.lt.0.d0 .or. xnu.ge.1.d0) call xerror ( 'd9knus xnu must 1 be ge 0 and lt 1', 33, 1, 2) if (x.le.0.) call xerror ( 'd9knus x must be gt 0', 22, 2, 2) c iswtch = 0 if (x.gt.2.0d0) go to 50 c c x is small. compute k-sub-xnu (x) and the derivative of k-sub-xnu (x) c then find k-sub-xnu+1 (x). xnu is reduced to the interval (-.5,+.5) c then to (0., .5), because k of negative order (-nu) = k of positive c order (+nu). c v = xnu if (xnu.gt.0.5d0) v = 1.0d0 - xnu c c carefully find (x/2)**xnu and z**xnu where z = x*x/4. alnz = 2.d0 * (dlog(x) - aln2) c if (x.gt.xnu) go to 20 if (-0.5d0*xnu*alnz-aln2-dlog(xnu).gt.alnbig) call xerror ( 'd9knu 1s x so small bessel k-sub-xnu overflows', 45, 3, 2) c 20 vlnz = v*alnz x2tov = dexp (0.5d0*vlnz) ztov = 0.0d0 if (vlnz.gt.alnsml) ztov = x2tov**2 c a0 = 0.5d0*dgamma(1.0d0+v) b0 = 0.5d0*dgamma(1.0d0-v) c0 = -euler if (ztov.gt.0.5d0 .and. v.gt.xnusml) c0 = -0.75d0 + 1 dcsevl ((8.0d0*v)*v-1.0d0, c0kcs, ntc0k) c if (ztov.le.0.5d0) alpha(1) = (a0-ztov*b0)/v if (ztov.gt.0.5d0) alpha(1) = c0 - alnz*(0.75d0 + 1 dcsevl (vlnz/0.35d0+1.0d0, znu1cs, ntznu1))*b0 beta(1) = -0.5d0*(a0+ztov*b0) c z = 0.0d0 if (x.gt.xsml) z = 0.25d0*x*x nterms = max1 (2.0, 11.0+(8.*sngl(alnz)-25.19-alneps) 1 /(4.28-sngl(alnz))) do 30 i=2,nterms xi = i - 1 a0 = a0/(xi*(xi-v)) b0 = b0/(xi*(xi+v)) alpha(i) = (alpha(i-1)+2.0d0*xi*a0)/(xi*(xi+v)) beta(i) = (xi-0.5d0*v)*alpha(i) - ztov*b0 30 continue c bknu = alpha(nterms) bknud = beta(nterms) do 40 ii=2,nterms i = nterms + 1 - ii bknu = alpha(i) + bknu*z bknud = beta(i) + bknud*z 40 continue c expx = dexp(x) bknu = expx*bknu/x2tov c if (-0.5d0*(xnu+1.d0)*alnz-2.0d0*aln2.gt.alnbig) iswtch = 1 if (iswtch.eq.1) return bknud = expx*bknud*2.0d0/(x2tov*x) c if (xnu.le.0.5d0) bknu1 = v*bknu/x - bknud if (xnu.le.0.5d0) return c bknu0 = bknu bknu = -v*bknu/x - bknud bknu1 = 2.0d0*xnu*bknu/x + bknu0 return c c x is large. find k-sub-xnu (x) and k-sub-xnu+1 (x) with y. l. luke-s c rational expansion. c 50 sqrtx = dsqrt(x) if (x.gt.1.0d0/xsml) go to 90 an = -0.60 - 1.02/sngl(x) bn = -0.27 - 0.53/sngl(x) nterms = min0 (32, max1 (3.0, an+bn*alneps)) c do 80 inu=1,2 xmu = 0.d0 if (inu.eq.1 .and. xnu.gt.xnusml) xmu = (4.0d0*xnu)*xnu if (inu.eq.2) xmu = 4.0d0*(dabs(xnu)+1.d0)**2 c a(1) = 1.0d0 - xmu a(2) = 9.0d0 - xmu a(3) = 25.0d0 - xmu if (a(2).eq.0.d0) result = sqpi2*(16.d0*x+xmu+7.d0) / 1 (16.d0*x*sqrtx) if (a(2).eq.0.d0) go to 70 c alpha(1) = 1.0d0 alpha(2) = (16.d0*x+a(2))/a(2) alpha(3) = ((768.d0*x+48.d0*a(3))*x + a(2)*a(3))/(a(2)*a(3)) c beta(1) = 1.0d0 beta(2) = (16.d0*x+(xmu+7.d0))/a(2) beta(3) = ((768.d0*x+48.d0*(xmu+23.d0))*x + 1 ((xmu+62.d0)*xmu+129.d0))/(a(2)*a(3)) c if (nterms.lt.4) go to 65 do 60 i=4,nterms n = i - 1 x2n = 2*n - 1 c a(i) = (x2n+2.d0)**2 - xmu qq = 16.d0*x2n/a(i) p1 = -x2n*(dble(float(12*n*n-20*n))-a(1))/((x2n-2.d0)*a(i)) 1 - qq*x p2 = (dble(float(12*n*n-28*n+8))-a(1))/a(i) - qq*x p3 = -x2n*a(i-3)/((x2n-2.d0)*a(i)) c alpha(i) = -p1*alpha(i-1) - p2*alpha(i-2) - p3*alpha(i-3) beta(i) = -p1*beta(i-1) - p2*beta(i-2) - p3*beta(i-3) 60 continue c 65 result = sqpi2*beta(nterms)/(sqrtx*alpha(nterms)) c 70 if (inu.eq.1) bknu = result if (inu.eq.2) bknu1 = result 80 continue return c 90 bknu = sqpi2/sqrtx bknu1 = bknu return c end 
Input File: f.omh/d9knus.f.omh
6.2: d9lgmc

6.2.a: Prototype
doublereal f2cad::d9lgmc_(doublereal *x);

6.2.b: Fortran Source
6.2.1: d9lgmc.f

6.2.c: Description
Change this paragraph to a description of this example / test.  # include <f2cad/d9lgmc.hpp> test_result d9lgmc(void) { // Replace the following line by an example / test of d9lgmc // that returns test_pass if the test passes and test_fail otherwise. return test_none; } 
Input File: example/d9lgmc.cpp
6.2.1: d9lgmc.f
 double precision function d9lgmc(x) c***begin prologue d9lgmc c***date written 770601 (yymmdd) c***revision date 820801 (yymmdd) c***category no. c7e c***keywords complete gamma function,correction factor, c double precision,gamma function,logarithm, c special function c***author fullerton, w., (lanl) c***purpose computes the d.p. log gamma correction factor for c x .ge. 10. so that dlog(dgamma(x)) = dlog(dsqrt(2*pi)) + c (x-5.)*dlog(x) - x + d9lgmc(x) c***description c c compute the log gamma correction factor for x .ge. 10. so that c dlog (dgamma(x)) = dlog(dsqrt(2*pi)) + (x-.5)*dlog(x) - x + d9lgmc(x) c c series for algm on the interval 0. to 1.00000e-02 c with weighted error 1.28e-31 c log weighted error 30.89 c significant figures required 29.81 c decimal places required 31.48 c***references (none) c***routines called d1mach,dcsevl,initds,xerror c***end prologue d9lgmc double precision x, algmcs(15), xbig, xmax, dcsevl, d1mach data algmcs( 1) / +.1666389480 4518632472 0572965082 2 d+0 / data algmcs( 2) / -.1384948176 0675638407 3298605913 5 d-4 / data algmcs( 3) / +.9810825646 9247294261 5717154748 7 d-8 / data algmcs( 4) / -.1809129475 5724941942 6330626671 9 d-10 / data algmcs( 5) / +.6221098041 8926052271 2601554341 6 d-13 / data algmcs( 6) / -.3399615005 4177219443 0333059966 6 d-15 / data algmcs( 7) / +.2683181998 4826987489 5753884666 6 d-17 / data algmcs( 8) / -.2868042435 3346432841 4462239999 9 d-19 / data algmcs( 9) / +.3962837061 0464348036 7930666666 6 d-21 / data algmcs( 10) / -.6831888753 9857668701 1199999999 9 d-23 / data algmcs( 11) / +.1429227355 9424981475 7333333333 3 d-24 / data algmcs( 12) / -.3547598158 1010705471 9999999999 9 d-26 / data algmcs( 13) / +.1025680058 0104709120 0000000000 0 d-27 / data algmcs( 14) / -.3401102254 3167487999 9999999999 9 d-29 / data algmcs( 15) / +.1276642195 6300629333 3333333333 3 d-30 / data nalgm, xbig, xmax / 0, 2*0.d0 / c***first executable statement d9lgmc if (nalgm.ne.0) go to 10 nalgm = initds (algmcs, 15, sngl(d1mach(3)) ) xbig = 1.0d0/dsqrt(d1mach(3)) xmax = dexp (dmin1(dlog(d1mach(2)/12.d0), -dlog(12.d0*d1mach(1)))) c 10 if (x.lt.10.d0) call xerror ( 'd9lgmc x must be ge 10', 23, 1, 2) if (x.ge.xmax) go to 20 c d9lgmc = 1.d0/(12.d0*x) if (x.lt.xbig) d9lgmc = dcsevl (2.0d0*(10.d0/x)**2-1.d0, algmcs, 1 nalgm) / x return c 20 d9lgmc = 0.d0 call xerror ( 'd9lgmc x so big d9lgmc underflows', 34, 2, 1) return c end 
Input File: f.omh/d9lgmc.f.omh
6.3: daxpy

6.3.a: Prototype
int f2cad::daxpy_(integer *n, doublereal *da, doublereal *dx, integer *incx, doublereal *dy, integer *incy);

6.3.b: Fortran Source
6.3.1: daxpy.f

6.3.c: Description
This example uses daxpy.f to compute   $\left( \begin{array}{c} f_0 \\ f_1 \end{array} \right) = 5 \left( \begin{array}{c} x_0 \\ x_1 \end{array} \right) + \left( \begin{array}{c} 3 \\ 4 \end{array} \right)$  Using the 8.3: f2cad_link routines f2cad::Independent and f2cad::Dependent, this defines the function   $f(x) = \left( \begin{array}{c} 5 x_0 + 3 \\ 5 x_1 + 4 \end{array} \right)$  We check that the derivative of this function, calculated using the f2cad::Partial routine, satisfies   $f^{(1)} (x) = \left( \begin{array}{cc} 5 & 0 \\ 0 & 5 \end{array} \right)$   # include <f2cad/daxpy.hpp> test_result daxpy(void) { bool ok = true; // Input values for daxpy doublereal a = 5.; integer n = 2; doublereal x[2]; x[0] = 1.; x[1] = 2.; integer incx = 1; integer m = 2; doublereal f[2]; f[0] = 3.; f[1] = 4.; integer incf = 1; // declare independent variables f2cad::Independent(n, x); // set f = a * x + f f2cad::daxpy_(&n, &a, x, &incx, f, &incf); // declare dependent variables f2cad::Dependent(m, f); double p; p = f2cad::Partial<doublereal>(0, 0); ok &= f2cad::near_equal(p, 5., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(1, 0); ok &= f2cad::near_equal(p, 0., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(0, 1); ok &= f2cad::near_equal(p, 0., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(1, 1); ok &= f2cad::near_equal(p, 5., 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/daxpy.cpp
6.3.1: daxpy.f
 subroutine daxpy(n,da,dx,incx,dy,incy) c***begin prologue daxpy c***date written 791001 (yymmdd) c***revision date 820801 (yymmdd) c***category no. d1a7 c***keywords blas,double precision,linear algebra,triad,vector c***author lawson, c. l., (jpl) c hanson, r. j., (snla) c kincaid, d. r., (u. of texas) c krogh, f. t., (jpl) c***purpose d.p computation y = a*x + y c***description c c b l a s subprogram c description of parameters c c --input-- c n number of elements in input vector(s) c da double precision scalar multiplier c dx double precision vector with n elements c incx storage spacing between elements of dx c dy double precision vector with n elements c incy storage spacing between elements of dy c c --output-- c dy double precision result (unchanged if n .le. 0) c c overwrite double precision dy with double precision da*dx + dy. c for i = 0 to n-1, replace dy(ly+i*incy) with da*dx(lx+i*incx) + c dy(ly+i*incy), where lx = 1 if incx .ge. 0, else lx = (-incx)*n c and ly is defined in a similar way using incy. c***references lawson c.l., hanson r.j., kincaid d.r., krogh f.t., c *basic linear algebra subprograms for fortran usage*, c algorithm no. 539, transactions on mathematical c software, volume 5, number 3, september 1979, 308-323 c***routines called (none) c***end prologue daxpy c double precision dx(1),dy(1),da c***first executable statement daxpy if(n.le.0.or.da.eq.0.d0) return if(incx.eq.incy) if(incx-1) 5,20,60 5 continue c c code for nonequal or nonpositive increments. c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop so remaining vector length is a multiple of 4. c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return c c code for equal, positive, nonunit increments. c 60 continue ns = n*incx do 70 i=1,ns,incx dy(i) = da*dx(i) + dy(i) 70 continue return end 
Input File: omh/daxpy.f.omh
6.4: dbesks

6.4.a: Prototype
int f2cad::dbesks_(doublereal *xnu, doublereal *x, integer *nin, doublereal *bk);

6.4.b: Fortran Source
6.4.1: dbesks.f

6.4.c: Notation
We use   K_{\nu} (x) to denote the modified Bessel function of the second kind; i.e., it is the solution of the ODE   $x^2 y^{(2)} (x) + x y^{(1)} (x) - ( x^2 + \nu^2 ) = 0$  that goes to zero as   x \rightarrow \infty .

6.4.d: Description
The routine dbesks computes the values   K_{\nu + i} (x) where   x > 0 ,   -1 < \nu < 1 ,   i = 0 , 1 , \ldots , nin-1 (if   nin > 0 ),   i = 0 , -1 , \ldots , nin+1 (if   nin < 0 ).

6.4.e: Special Case
By Formula 10.2.17 of Abramowitz and Stegun   $K_{3/2} ( x ) = \sqrt{ \pi / ( 2 x ) } ( 1 + 1 / x ) \exp( - x )$   # include <cmath> # include <f2cad/dbesks.hpp> test_result dbesks(void) { bool ok = true; // independent variable integer n = 1; // number of independent variables doublereal x[1]; // vector of independent variables x[0] = .75; // value of independent variable f2cad::Independent(n, x); // declare x as independent variable vector // call the routine doublereal xnu = .5; doublereal bk[2]; integer nin = 2; f2cad::dbesks_(&xnu, x+0, &nin, bk); // dependent variables integer m = 1; // number of dependent variables doublereal f[1]; // vector of dependent variables f[0] = bk[1]; // value f[0] = K_(nu+1) (x[0]) f2cad::Dependent(m, f); // declare f as dependent variable vector // check function value double x0 = f2cad::Value( x[0] ); double f0 = f2cad::Value( f[0] ); double pi = 4. * std::atan(1.); double term = 1. / sqrt( 2. * x0 / pi ); double check = term * (1. + 1. / x0 ) * exp(-x0); ok &= f2cad::near_equal(f0, check, 1e-10, 1e-10); // Evaluate the derivative of f[0] w.r.t x[0] double p = f2cad::Partial<doublereal>(0, 0); check = - check; check += - term * exp(-x0) / (x0 * x0); check += - (1. + 1. / x0) * exp(-x0) * term * term * term / pi; ok &= f2cad::near_equal(p, check, 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/dbesks.cpp
6.4.1: dbesks.f
 subroutine dbesks(xnu,x,nin,bk) c***begin prologue dbesks c***date written 770601 (yymmdd) c***revision date 820801 (yymmdd) c***category no. c10b3 c***keywords bessel function,double precision,fractional order, c modified bessel function,sequence,special function, c third kind c***author fullerton, w., (lanl) c***purpose computes a d.p. sequence of modified bessel functions of c the third kind of fractional order. c***description c c dbesks computes a sequence of modified bessel functions of the third c kind of order xnu + i at x, where x .gt. 0, xnu lies in (-1,1), c and i = 0, 1, ... , nin - 1, if nin is positive and i = 0, 1, ... , c nin + 1, if nin is negative. on return, the vector bk(.) contains c the results at x for order starting at xnu. xnu, x, and bk are c double precision. nin is an integer. c***references (none) c***routines called d1mach,dbskes,xerror c***end prologue dbesks double precision xnu, x, bk(1), expxi, xmax, d1mach data xmax / 0.d0 / c***first executable statement dbesks if (xmax.eq.0.d0) xmax = -dlog (d1mach(1)) c if (x.gt.xmax) call xerror ( 'dbesks x so big bessel k underflows 1', 36, 1, 2) c call dbskes (xnu, x, nin, bk) c expxi = dexp (-x) n = iabs (nin) do 20 i=1,n bk(i) = expxi * bk(i) 20 continue c return end 
Input File: f.omh/dbesks.f.omh
6.5: dbskes

6.5.a: Prototype
int f2cad::dbskes_(doublereal *xnu, doublereal *x, integer *nin, doublereal *bke);

6.5.b: Fortran Source
6.5.1: dbskes.f

6.5.c: Notation
We use   K_{\nu} (x) to denote the modified Bessel function of the second kind; i.e., it is the solution of the ODE   $x^2 y^{(2)} (x) + x y^{(1)} (x) - ( x^2 + \nu^2 ) = 0$  that goes to zero as   x \rightarrow \infty .

6.5.d: Description
The routine dbskes computes the values   \exp(x) K_{\nu + i} (x) where   x > 0 ,   -1 < \nu < 1 ,   i = 0 , 1 , \ldots , nin-1 (if   nin > 0 ),   i = 0 , -1 , \ldots , nin+1 (if   nin < 0 ).

6.5.e: Special Case
By Formula 10.2.17 of Abramowitz and Stegun   $K_{3/2} ( x ) = \sqrt{ \pi / ( 2 x ) } ( 1 + 1 / x ) \exp( - x )$   # include <cmath> # include <f2cad/dbskes.hpp> test_result dbskes(void) { bool ok = true; // independent variable integer n = 1; // number of independent variables doublereal x[1]; // vector of independent variables x[0] = .75; // value of independent variable f2cad::Independent(n, x); // declare x as independent variable vector // call the routine doublereal xnu = .5; doublereal bke[2]; integer nin = 2; f2cad::dbskes_(&xnu, x+0, &nin, bke); // dependent variables integer m = 1; // number of dependent variables doublereal f[1]; // vector of dependent variables f[0] = bke[1]; // value f[0] = K_(nu+1) (x[0]) * exp( x[0] ) f2cad::Dependent(m, f); // declare f as dependent variable vector // check function value double x0 = f2cad::Value( x[0] ); double f0 = f2cad::Value( f[0] ); double pi = 4. * std::atan(1.); double term = 1. / sqrt( 2. * x0 / pi ); double check = term * (1. + 1. / x0 ); ok &= f2cad::near_equal(f0, check, 1e-10, 1e-10); // Evaluate the derivative of f[0] w.r.t x[0] double p = f2cad::Partial<doublereal>(0, 0); check = - (1. + 1. / x0) * term * term * term / pi; check += - term / (x0 * x0); ok &= f2cad::near_equal(p, check, 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/dbskes.cpp
6.5.1: dbskes.f
 subroutine dbskes(xnu,x,nin,bke) c***begin prologue dbskes c***date written 770601 (yymmdd) c***revision date 820801 (yymmdd) c***category no. c10b3 c***keywords bessel function,double precision,exponentially scaled, c fractional order,modified bessel function,sequence, c special function,third kind c***author fullerton, w., (lanl) c***purpose computes a d.p. sequence of exponentially scaled modified c bessel functions of the third kind of fractional order. c***description c c dbskes(xnu,x,nin,bke) computes a double precision sequence c of exponentially scaled modified bessel functions c of the third kind of order xnu + i at x, where x .gt. 0, c xnu lies in (-1,1), and i = 0, 1, ... , nin - 1, if nin is positive c and i = 0, -1, ... , nin + 1, if nin is negative. on return, the c vector bke(.) contains the results at x for order starting at xnu. c xnu, x, and bke are double precison. nin is integer. c***references (none) c***routines called d1mach,d9knus,xerror c***end prologue dbskes double precision xnu, x, bke(1), bknu1, v, vincr, vend, alnbig, 1 d1mach data alnbig / 0.d0 / c***first executable statement dbskes if (alnbig.eq.0.d0) alnbig = dlog (d1mach(2)) c v = dabs(xnu) n = iabs(nin) c if (v.ge.1.d0) call xerror ( 'dbskes abs(xnu) must be lt 1', 29, 1 2, 2) if (x.le.0.d0) call xerror ( 'dbskes x is le 0', 17, 3, 2) if (n.eq.0) call xerror ( 'dbskes n the number in the sequence is 1 0', 41, 4, 2) c call d9knus (v, x, bke(1), bknu1, iswtch) if (n.eq.1) return c vincr = sign (1.0, float(nin)) direct = vincr if (xnu.ne.0.d0) direct = vincr*dsign(1.d0, xnu) if (iswtch.eq.1 .and. direct.gt.0.) call xerror ( 'dbskes x so sm 1all bessel k-sub-xnu+1 overflows', 47, 5, 2) bke(2) = bknu1 c if (direct.lt.0.) call d9knus (dabs(xnu+vincr), x, bke(2), bknu1, 1 iswtch) if (n.eq.2) return c vend = dabs (xnu+dble(float(nin))) - 1.0d0 if ((vend-.5d0)*dlog(vend)+0.27d0-vend*(dlog(x)-.694d0).gt.alnbig) 1 call xerror ( 'dbskes x so small or abs(nu) so big that bessel 2 k-sub-nu overflows', 67, 5, 2) c v = xnu do 10 i=3,n v = v + vincr bke(i) = 2.0d0*v*bke(i-1)/x + bke(i-2) 10 continue c return end 
Input File: f.omh/dbskes.f.omh
6.6: dcsevl

6.6.a: Prototype
doublereal f2cad::dcsevl_(doublereal *x, doublereal *a, integer *n);

6.6.b: Fortran Source
6.6.1: dcsevl.f

6.6.c: Description
The routine dcsevl computes values for a Chebyshev Series; i.e., the return value is   $f(x) = a_0 / 2 + \sum_{i=1}^{n-1} a_i T_i (x)$  where   T_i (x) is the i-th order Chebyshev polynomial of the first kind; i.e., it is defined by the recurrence relation   $\begin{array}{rcl} T_0 (x) & = & 1 \\ T_1 (x) & = & x \\ T_{i+1} (x) & = & 2 x T_i (x) - T_{i-1} (x) \end{array}$  It follows that the derivatives   T_i^{(1)} (x) satisfy the recurrence relation   $\begin{array}{rcl} T_0^{(1)} (x) & = & 0 \\ T_1^{(1)} (x) & = & 1 \\ T_{i+1}^{(1)} (x) & = & 2 T_i (x) + 2 x T_i^{(1)} (x) - T_{i-1}^{(1)} (x) \end{array}$  and the derivative of   f(x) is given by   $f^{(1)} (x) = \sum_{i=1}^{n-1} a_i T_i^{(1)} (x)$   # include <f2cad/dcsevl.hpp> test_result dcsevl(void) { bool ok = true; // independent variable integer N = 1; double x = .5; doublereal X[1]; X[0] = x; // dependent variable integer M = 1; doublereal F[1]; // polynomial coefficients integer n = 3; double a[3]; doublereal A[3]; for(integer i = 0; i < n; i++) { a[i] = 1. / double(i+1); A[i] = a[i]; } // declare independent variables f2cad::Independent(N, X); F[0] = f2cad::dcsevl_(X, A, &n); // declare dependent variables f2cad::Dependent(M, F); // evaluate the Chebyshev polynomials double t_0 = 1.; double t_1 = x; double t_2 = 2. * x * t_1 - t_0; // check function value double value = f2cad::Value( F[0] ); double check = a[0] * t_0 / 2. + a[1] * t_1 + a[2] * t_2; ok &= f2cad::near_equal(check, value, 1e-10, 1e-10); // evaluate the Derivatives double dt_0 = 0.; double dt_1 = 1.; double dt_2 = 2. * t_1 + 2. * x * dt_1 - dt_0; value = f2cad::Partial<doublereal>(0, 0); check = a[1] * dt_1 + a[2] * dt_2; ok &= f2cad::near_equal(check, value, 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/dcsevl.cpp
6.6.1: dcsevl.f
 double precision function dcsevl(x,a,n) c***begin prologue dcsevl c***date written 770401 (yymmdd) c***revision date 820801 (yymmdd) c***category no. c3a2 c***keywords chebyshev,fnlib,special function c***author fullerton, w., (lanl) c***purpose evaluate the double precision n-term chebyshev series a c at x. c***description c c evaluate the n-term chebyshev series a at x. adapted from c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). c w. fullerton, c-3, los alamos scientific laboratory. c c input arguments -- c x double precision value at which the series is to be evaluated. c a double precision array of n terms of a chebyshev series. in c evaluating a, only half of the first coefficient is summed. c n number of terms in array a. c***references (none) c***routines called xerror c***end prologue dcsevl c double precision a(n),x,twox,b0,b1,b2 c***first executable statement dcsevl if(n.lt.1)call xerror( 'dcsevl number of terms le 0', 28, 2,2) if(n.gt.1000) call xerror ( 'dcsevl number of terms gt 1000', 1 31, 3, 2) if ((x.lt.-1.d0) .or. (x.gt.1.d0)) call xerror ( 'dcsevl x outsi 1de (-1,+1)', 25, 1, 1) c twox = 2.0d0*x b1 = 0.d0 b0=0.d0 do 10 i=1,n b2=b1 b1=b0 ni = n - i + 1 b0 = twox*b1 - b2 + a(ni) 10 continue c dcsevl = 0.5d0 * (b0-b2) c return end 
Input File: omh/dcsevl.f.omh
6.7: ddot

6.7.a: Prototype
doublereal f2cad::ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy);

6.7.b: Fortran Source
6.7.1: ddot.f

6.7.c: Description
This example uses the routine ddot.f to compute   $f_0 = \left( \begin{array}{cc} x_0 & x_1 \end{array} \right) \left( \begin{array}{c} x_2 \\ x_3 \end{array} \right)$  Using the 8.3: f2cad_link routines f2cad::Independent and f2cad::Dependent, this defines the function   $f(x) = x_0 * x_2 + x_1 * x_3$  We check that the derivative of this function, calculated using the f2cad::Partial routine, satisfies   $f^{(1)} (x) = \left( \begin{array}{cccc} x_2 & x_3 & x_0 & x_1 \end{array} \right)$   # include <f2cad/ddot.hpp> test_result ddot(void) { bool ok = true; integer n = 4; doublereal x[4]; x[0] = 1.; x[1] = 2.; x[2] = 3.; x[3] = 4.; // declare independent variables f2cad::Independent(n, x); // set f = x[0] * x[2] + x[1] * x[3] integer m = 1; doublereal f[1]; integer incx = 1; integer n2 = n / 2; f[0] = f2cad::ddot_ (&n2, x, &incx, x+n2, &incx); // declare dependent variables f2cad::Dependent(m, f); double p; p = f2cad::Partial<doublereal>(0, 0); ok &= f2cad::near_equal(p, f2cad::Value(x[2]), 1e-10, 1e-10); p = f2cad::Partial<doublereal>(0, 1); ok &= f2cad::near_equal(p, f2cad::Value(x[3]), 1e-10, 1e-10); p = f2cad::Partial<doublereal>(0, 2); ok &= f2cad::near_equal(p, f2cad::Value(x[0]), 1e-10, 1e-10); p = f2cad::Partial<doublereal>(0, 3); ok &= f2cad::near_equal(p, f2cad::Value(x[1]), 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/ddot.cpp
6.7.1: ddot.f
 double precision function ddot(n,dx,incx,dy,incy) c***begin prologue ddot c***date written 791001 (yymmdd) c***revision date 820801 (yymmdd) c***category no. d1a4 c***keywords blas,double precision,inner product,linear algebra,vector c***author lawson, c. l., (jpl) c hanson, r. j., (snla) c kincaid, d. r., (u. of texas) c krogh, f. t., (jpl) c***purpose d.p. inner product of d.p. vectors c***description c c b l a s subprogram c description of parameters c c --input-- c n number of elements in input vector(s) c dx double precision vector with n elements c incx storage spacing between elements of dx c dy double precision vector with n elements c incy storage spacing between elements of dy c c --output-- c ddot double precision dot product (zero if n .le. 0) c c returns the dot product of double precision dx and dy. c ddot = sum for i = 0 to n-1 of dx(lx+i*incx) * dy(ly+i*incy) c where lx = 1 if incx .ge. 0, else lx = (-incx)*n, and ly is c defined in a similar way using incy. c***references lawson c.l., hanson r.j., kincaid d.r., krogh f.t., c *basic linear algebra subprograms for fortran usage*, c algorithm no. 539, transactions on mathematical c software, volume 5, number 3, september 1979, 308-323 c***routines called (none) c***end prologue ddot c double precision dx(1),dy(1) c***first executable statement ddot ddot = 0.d0 if(n.le.0)return if(incx.eq.incy) if(incx-1) 5,20,60 5 continue c c code for unequal or nonpositive increments. c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n ddot = ddot + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1. c c c clean-up loop so remaining vector length is a multiple of 5. c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m ddot = ddot + dx(i)*dy(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 ddot = ddot + dx(i)*dy(i) + dx(i+1)*dy(i+1) + 1 dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue return c c code for positive equal increments .ne.1. c 60 continue ns = n*incx do 70 i=1,ns,incx ddot = ddot + dx(i)*dy(i) 70 continue return end 
Input File: omh/ddot.f.omh
6.8: dgamlm

6.8.a: Prototype
int f2cad::dgamlm_(doublereal *xmin, doublereal *xmax);

6.8.b: Fortran Source
6.8.1: dgamlm.f

6.8.c: Description
This routine computes the allowable limits for the argument to the Gamma function. These limits should also work for the exponential function.  # include <f2cad/dgamlm.hpp> test_result dgamlm(void) { bool ok = true; doublereal ad_xmin, ad_xmax; f2cad::dgamlm_(&ad_xmin, &ad_xmax); double xmin = f2cad::Value(ad_xmin); double xmax = f2cad::Value(ad_xmax); ok &= xmin < 0.; ok &= xmax > 0.; ok &= (xmax - xmin) > 50.; ok &= f2cad::near_equal( log(exp(xmin)), xmin, 1e-10, 1e-10); ok &= f2cad::near_equal( log(exp(xmax)), xmax, 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/dgamlm.cpp
6.8.1: dgamlm.f
 subroutine dgamlm(xmin,xmax) c***begin prologue dgamlm c***date written 770601 (yymmdd) c***revision date 820801 (yymmdd) c***category no. c7a,r2 c***keywords complete gamma function,double precision,gamma function, c limits,special function c***author fullerton, w., (lanl) c***purpose computes the d.p. minimum and maximum bounds for x in c gamma(x). c***description c c calculate the minimum and maximum legal bounds for x in gamma(x). c xmin and xmax are not the only bounds, but they are the only non- c trivial ones to calculate. c c output arguments -- c xmin double precision minimum legal value of x in gamma(x). any c smaller value of x might result in underflow. c xmax double precision maximum legal value of x in gamma(x). any c larger value of x might cause overflow. c***references (none) c***routines called d1mach,xerror c***end prologue dgamlm double precision xmin, xmax, alnbig, alnsml, xln, xold, d1mach c***first executable statement dgamlm alnsml = dlog(d1mach(1)) xmin = -alnsml do 10 i=1,10 xold = xmin xln = dlog(xmin) xmin = xmin - xmin*((xmin+0.5d0)*xln - xmin - 0.2258d0 + alnsml) 1 / (xmin*xln+0.5d0) if (dabs(xmin-xold).lt.0.005d0) go to 20 10 continue call xerror ( 'dgamlm unable to find xmin', 27, 1, 2) c 20 xmin = -xmin + 0.01d0 c alnbig = dlog (d1mach(2)) xmax = alnbig do 30 i=1,10 xold = xmax xln = dlog(xmax) xmax = xmax - xmax*((xmax-0.5d0)*xln - xmax + 0.9189d0 - alnbig) 1 / (xmax*xln-0.5d0) if (dabs(xmax-xold).lt.0.005d0) go to 40 30 continue call xerror ( 'dgamlm unable to find xmax', 27, 2, 2) c 40 xmax = xmax - 0.01d0 xmin = dmax1 (xmin, -xmax+1.d0) c return end 
Input File: omh/dgamlm.f.omh
6.9: dgamma

6.9.a: Prototype
doublereal f2cad::dgamma_(doublereal *x);

6.9.b: Fortran Source
6.9.1: dgamma.f

6.9.c: Description
For arbitrary   x and positive integer values of   \ell ,   $\begin{array}{rcl} \Gamma(x) & = & \int_0^\infty t^{x - 1} \exp( - t ) \; {\bf d} t \\ \Gamma(\ell+1) & = & \ell ! \\ \Gamma^{(1)} (\ell+1) & = & \ell ! \left( - \gamma + \sum_{k=1}^\ell \frac{1}{k} \right) \end{array}$  where   \gamma is the Euler-Mascheroni constant   $\gamma = 0.57721 56649 01532 86060 65120 90082 ...$   # include <f2cad/dgamma.hpp> test_result dgamma(void) { bool ok = true; // independent variables integer n = 1; // number of independent variables double ell = 3.; doublereal x[1]; // vector of independent variables x[0] = ell+1.; // value of independent variables f2cad::Independent(n, x); // declare x as independent variable vector // dependent variables integer m = 1; // number of dependent variables doublereal f[1]; // vector of dependent variables f[0] = f2cad::dgamma_(x+0); // compute value f[0] = Gamma(x[0]) f2cad::Dependent(m, f); // declare f as dependent variable vector // check the value of f ok &= f2cad::near_equal(f2cad::Value(f[0]), 6., 1e-10, 1e-10); // Compute derivative of f w.r.t x using f2cad common interface. // This general interface works with adolc, cppad, and fadbad. // See GetStarted for package specific more efficient code. double gamma = 0.577215664901532; double p = f2cad::Partial<doublereal>(0, 0); double check = 6. * ( - gamma + 1. + 1./2. + 1./3. ); ok &= f2cad::near_equal(p, check, 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/dgamma.cpp
6.9.1: dgamma.f
 double precision function dgamma(x) c***begin prologue dgamma c***date written 770601 (yymmdd) c***revision date 820801 (yymmdd) c***category no. c7a c***keywords complete gamma function,double precision,gamma function, c special function c***author fullerton, w., (lanl) c***purpose computes the d.p. complete gamma function. c***description c c dgamma(x) calculates the double precision complete gamma function c for double precision argument x. c c series for gam on the interval 0. to 1.00000e+00 c with weighted error 5.79e-32 c log weighted error 31.24 c significant figures required 30.00 c decimal places required 32.05 c***references (none) c***routines called d1mach,d9lgmc,dcsevl,dgamlm,dint,initds,xerror c***end prologue dgamma double precision x, gamcs(42), dxrel, pi, sinpiy, sq2pil, xmax, 1 xmin, y, dint, d9lgmc, dcsevl, d1mach c data gam cs( 1) / +.8571195590 9893314219 2006239994 2 d-2 / data gam cs( 2) / +.4415381324 8410067571 9131577165 2 d-2 / data gam cs( 3) / +.5685043681 5993633786 3266458878 9 d-1 / data gam cs( 4) / -.4219835396 4185605010 1250018662 4 d-2 / data gam cs( 5) / +.1326808181 2124602205 8400679635 2 d-2 / data gam cs( 6) / -.1893024529 7988804325 2394702388 6 d-3 / data gam cs( 7) / +.3606925327 4412452565 7808221722 5 d-4 / data gam cs( 8) / -.6056761904 4608642184 8554829036 5 d-5 / data gam cs( 9) / +.1055829546 3022833447 3182350909 3 d-5 / data gam cs( 10) / -.1811967365 5423840482 9185589116 6 d-6 / data gam cs( 11) / +.3117724964 7153222777 9025459316 9 d-7 / data gam cs( 12) / -.5354219639 0196871408 7408102434 7 d-8 / data gam cs( 13) / +.9193275519 8595889468 8778682594 0 d-9 / data gam cs( 14) / -.1577941280 2883397617 6742327395 3 d-9 / data gam cs( 15) / +.2707980622 9349545432 6654043308 9 d-10 / data gam cs( 16) / -.4646818653 8257301440 8166105893 3 d-11 / data gam cs( 17) / +.7973350192 0074196564 6076717535 9 d-12 / data gam cs( 18) / -.1368078209 8309160257 9949917230 9 d-12 / data gam cs( 19) / +.2347319486 5638006572 3347177168 8 d-13 / data gam cs( 20) / -.4027432614 9490669327 6657053469 9 d-14 / data gam cs( 21) / +.6910051747 3721009121 3833697525 7 d-15 / data gam cs( 22) / -.1185584500 2219929070 5238712619 2 d-15 / data gam cs( 23) / +.2034148542 4963739552 0102605193 2 d-16 / data gam cs( 24) / -.3490054341 7174058492 7401294910 8 d-17 / data gam cs( 25) / +.5987993856 4853055671 3505106602 6 d-18 / data gam cs( 26) / -.1027378057 8722280744 9006977843 1 d-18 / data gam cs( 27) / +.1762702816 0605298249 4275966074 8 d-19 / data gam cs( 28) / -.3024320653 7353062609 5877211204 2 d-20 / data gam cs( 29) / +.5188914660 2183978397 1783355050 6 d-21 / data gam cs( 30) / -.8902770842 4565766924 4925160106 6 d-22 / data gam cs( 31) / +.1527474068 4933426022 7459689130 6 d-22 / data gam cs( 32) / -.2620731256 1873629002 5732833279 9 d-23 / data gam cs( 33) / +.4496464047 8305386703 3104657066 6 d-24 / data gam cs( 34) / -.7714712731 3368779117 0390152533 3 d-25 / data gam cs( 35) / +.1323635453 1260440364 8657271466 6 d-25 / data gam cs( 36) / -.2270999412 9429288167 0231381333 3 d-26 / data gam cs( 37) / +.3896418998 0039914493 2081663999 9 d-27 / data gam cs( 38) / -.6685198115 1259533277 9212799999 9 d-28 / data gam cs( 39) / +.1146998663 1400243843 4761386666 6 d-28 / data gam cs( 40) / -.1967938586 3451346772 9510399999 9 d-29 / data gam cs( 41) / +.3376448816 5853380903 3489066666 6 d-30 / data gam cs( 42) / -.5793070335 7821357846 2549333333 3 d-31 / data pi / 3.1415926535 8979323846 2643383279 50 d0 / data sq2pil / 0.9189385332 0467274178 0329736405 62 d0 / data ngam, xmin, xmax, dxrel / 0, 3*0.d0 / c***first executable statement dgamma if (ngam.ne.0) go to 10 ngam = initds (gamcs, 42, 0.1*sngl(d1mach(3)) ) c call dgamlm (xmin, xmax) dxrel = dsqrt (d1mach(4)) c 10 y = dabs(x) if (y.gt.10.d0) go to 50 c c compute gamma(x) for -xbnd .le. x .le. xbnd. reduce interval and find c gamma(1+y) for 0.0 .le. y .lt. 1.0 first of all. c n = x if (x.lt.0.d0) n = n - 1 y = x - dble(float(n)) n = n - 1 dgamma = 0.9375d0 + dcsevl (2.d0*y-1.d0, gamcs, ngam) if (n.eq.0) return c if (n.gt.0) go to 30 c c compute gamma(x) for x .lt. 1.0 c n = -n if (x.eq.0.d0) call xerror ( 'dgamma x is 0', 14, 4, 2) if (x.lt.0.0 .and. x+dble(float(n-2)).eq.0.d0) call xerror ( 'dgam 1ma x is a negative integer', 31, 4, 2) if (x.lt.(-0.5d0) .and. dabs((x-dint(x-0.5d0))/x).lt.dxrel) call 1 xerror ( 'dgamma answer lt half precision because x too near ne 2gative integer', 68, 1, 1) c do 20 i=1,n dgamma = dgamma/(x+dble(float(i-1)) ) 20 continue return c c gamma(x) for x .ge. 2.0 and x .le. 10.0 c 30 do 40 i=1,n dgamma = (y+dble(float(i))) * dgamma 40 continue return c c gamma(x) for dabs(x) .gt. 10.0. recall y = dabs(x). c 50 if (x.gt.xmax) call xerror ( 'dgamma x so big gamma overflows', 1 32, 3, 2) c dgamma = 0.d0 if (x.lt.xmin) call xerror ( 'dgamma x so small gamma underflows' 1 , 35, 2, 1) if (x.lt.xmin) return c dgamma = dexp ((y-0.5d0)*dlog(y) - y + sq2pil + d9lgmc(y) ) if (x.gt.0.d0) return c if (dabs((x-dint(x-0.5d0))/x).lt.dxrel) call xerror ( 'dgamma ans 1wer lt half precision, x too near negative integer' , 61, 1, 1) c sinpiy = dsin (pi*y) if (sinpiy.eq.0.d0) call xerror ( 'dgamma x is a negative integer 1', 31, 4, 2) c dgamma = -pi/(y*sinpiy*dgamma) c return end 
Input File: f.omh/dgamma.f.omh
6.10: dgefa

6.10.a: Prototype
int f2cad::dgefa_(doublereal *a, integer *lda, integer *n, integer *ipvt, integer *info); int f2cad::dgesl_(doublereal *a, integer *lda, integer *n, integer *ipvt, doublereal *b, integer *job);

6.10.b: Fortran Source
6.10.1: dgefa.f

6.10.c: Description
This example uses the Linpack routines dgefa.f and dgesl.f to solve the linear equation   $\left( \begin{array}{cc} 1 & 2 \\ 0 & 3 \end{array} \right) \left( \begin{array}{c} f_0 \\ f_1 \end{array} \right) = \left( \begin{array}{c} b_0 \\ b_1 \end{array} \right)$  Using the 8.3: f2cad_link routines f2cad::Independent and f2cad::Dependent, this defines the function   $f(b) = \left( \begin{array}{c} b_0 - 2 b_1 / 3 \\ b_1 / 3 \end{array} \right)$  We check that the derivative of this function, calculated using the f2cad::Partial routine, satisfies   $f^{(1)} (b) = \left( \begin{array}{cc} 1 & - 2 / 3 \\ 0 & 1 / 3 \end{array} \right)$   # include <f2cad/dgefa.hpp> # include <f2cad/dgesl.hpp> test_result dgefa(void) { bool ok = true; // A is a 2 by 2 matrix in column major order doublereal A[4]; A[0] = doublereal(1); A[2] = doublereal(2); A[1] = doublereal(0); A[3] = doublereal(3); // b is a column vector of length 2 doublereal b[2]; b[0] = doublereal(1); b[1] = doublereal(2); // declare independent variables integer n = 2; f2cad::Independent(n, b); // f is a column vector of length 2 doublereal f[2]; integer i; for(i = 0; i < n; i++) f[i] = b[i]; // solve A * f = b integer ipvt[2]; integer lda = n; integer job = 0; integer info; f2cad::dgefa_ (A, &lda, &n, ipvt, &info); f2cad::dgesl_ (A, &lda, &n, ipvt, f, &job); // check info flag ok &= (info == 0); // construct the AD function object F : b -> f // f[1] = b[1] / 3 // f[0] = b[0] - 2 * b[1] / 3 integer m = 2; f2cad::Dependent(m, f); double p; p = f2cad::Partial<doublereal>(0, 0); ok &= f2cad::near_equal(p, 1., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(0, 1); ok &= f2cad::near_equal(p, -2./3., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(1, 0); ok &= f2cad::near_equal(p, 0., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(1, 1); ok &= f2cad::near_equal(p, 1./3., 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/dgefa.cpp
6.10.1: dgefa.f
 subroutine dgefa(a,lda,n,ipvt,info) c***begin prologue dgefa c***date written 780814 (yymmdd) c***revision date 820801 (yymmdd) c***category no. d2a1 c***keywords double precision,factor,linear algebra,linpack,matrix c***author moler, c. b., (u. of new mexico) c***purpose factors a double precision matrix by gaussian elimination. c***description c c dgefa factors a double precision matrix by gaussian elimination. c c dgefa is usually called by dgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for dgeco) = (1 + 9/n)*(time for dgefa) . c c on entry c c a double precision(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that dgesl or dgedi will divide by zero c if called. use rcond in dgeco for a reliable c indication of singularity. c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas daxpy,dscal,idamax c***references dongarra j.j., bunch j.r., moler c.b., stewart g.w., c *linpack users guide*, siam, 1979. c***routines called daxpy,dscal,idamax c***end prologue dgefa integer lda,n,ipvt(1),info double precision a(lda,1) c double precision t integer idamax,j,k,kp1,l,nm1 c c gaussian elimination with partial pivoting c c***first executable statement dgefa info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = idamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0.0d0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -1.0d0/a(k,k) call dscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0.0d0) info = n return end 
Input File: omh/dgefa.f.omh
6.11: dgesl

6.11.a: Prototype
int f2cad::dgefa_(doublereal *a, integer *lda, integer *n, integer *ipvt, integer *info); int f2cad::dgesl_(doublereal *a, integer *lda, integer *n, integer *ipvt, doublereal *b, integer *job);

6.11.b: Fortran Source
6.11.1: dgesl.f

6.11.c: Description
This example uses the Linpack routines dgefa.f and dgesl.f to solve the linear equation   $\left( \begin{array}{cc} 1 & 2 \\ 3 & 4 \end{array} \right) \left( \begin{array}{c} f_0 \\ f_1 \end{array} \right) = \left( \begin{array}{c} b_0 \\ b_1 \end{array} \right)$  Multiplying the top row by 3 and subtracting it from the bottom row we obtain the equivalent equation:   $\left( \begin{array}{cc} 1 & 2 \\ 0 & -2 \end{array} \right) \left( \begin{array}{c} f_0 \\ f_1 \end{array} \right) = \left( \begin{array}{c} b_0 \\ b_1 - 3 b_0 \end{array} \right)$  It follows that   $\left( \begin{array}{c} f_0 \\ f_1 \end{array} \right) = \left( \begin{array}{c} b_1 - 2 b_0 \\ 3 b_0 / 2 - b_1/ 2 \end{array} \right)$  The code below uses the 8.3: f2cad_link routines f2cad::Independent and f2cad::Dependent, to define the function   $f(b) = \left( \begin{array}{c} b_1 - 2 b_0 \\ 3 b_0 / 2 - b_1/ 2 \end{array} \right)$  We check that the derivative of this function, calculated using the f2cad::Partial routine, satisfies   $f^{(1)} (b) = \left( \begin{array}{cc} -2 & 1 \\ 3/2 & -1 / 2 \end{array} \right)$   # include <f2cad/dgefa.hpp> # include <f2cad/dgesl.hpp> test_result dgesl(void) { bool ok = true; // A is a 2 by 2 matrix in column major order doublereal A[4]; A[0] = doublereal(1); A[2] = doublereal(2); A[1] = doublereal(3); A[3] = doublereal(4); // b is a column vector of length 2 doublereal b[2]; b[0] = doublereal(1); b[1] = doublereal(2); // declare independent variables integer n = 2; f2cad::Independent(n, b); // f is a column vector of length 2 doublereal f[2]; integer i; for(i = 0; i < n; i++) f[i] = b[i]; // solve A * f = b integer ipvt[2]; integer lda = n; integer job = 0; integer info; f2cad::dgefa_ (A, &lda, &n, ipvt, &info); f2cad::dgesl_ (A, &lda, &n, ipvt, f, &job); // check info flag ok &= (info == 0); // construct the AD function object F : b -> f // f[0] = b[1] - 2 b[0] // f[1] = 3 b[0] / 2 - b[1]/ 2 integer m = 2; f2cad::Dependent(m, f); double p; p = f2cad::Partial<doublereal>(0, 0); ok &= f2cad::near_equal(p, -2., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(0, 1); ok &= f2cad::near_equal(p, 1., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(1, 0); ok &= f2cad::near_equal(p, 3./2., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(1, 1); ok &= f2cad::near_equal(p, -1./2., 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/dgesl.cpp
6.11.1: dgesl.f
 subroutine dgesl(a,lda,n,ipvt,b,job) c***begin prologue dgesl c***date written 780814 (yymmdd) c***revision date 820801 (yymmdd) c***category no. d2a1 c***keywords double precision,linear algebra,linpack,matrix,solve c***author moler, c. b., (u. of new mexico) c***purpose solves the double precision system a*x=b or trans(a)*x=b c using the factors computed by dgeco or dgefa. c***description c c dgesl solves the double precision system c a * x = b or trans(a) * x = b c using the factors computed by dgeco or dgefa. c c on entry c c a double precision(lda, n) c the output from dgeco or dgefa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c ipvt integer(n) c the pivot vector from dgeco or dgefa. c c b double precision(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if dgeco has set rcond .gt. 0.0 c or dgefa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call dgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call dgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas daxpy,ddot c***references dongarra j.j., bunch j.r., moler c.b., stewart g.w., c *linpack users guide*, siam, 1979. c***routines called daxpy,ddot c***end prologue dgesl integer lda,n,ipvt(1),job double precision a(lda,1),b(1) c double precision ddot,t integer k,kb,l,nm1 c***first executable statement dgesl nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call daxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call daxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n t = ddot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue c c now solve trans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end 
Input File: omh/dgesl.f.omh
6.12: dint

6.12.a: Prototype
doublereal f2cad::dint_(doublereal *x);

6.12.b: Fortran Source
6.12.1: dint.f

6.12.c: Description
This converts from a doublereal to an integer by rounding toward zero.  # include <f2cad/dint.hpp> test_result dint(void) { bool ok = true; // independent variables integer n = 2; // number of independent variables doublereal x[2]; // vector of independent variables x[0] = 4.5; // value of independent variables x[1] = -3.5; f2cad::Independent(n, x); // declare x as independent variable vector // dependent variables integer m = n; // number of dependent variables doublereal f[2]; // vector of dependent variables integer i; for(i = 0; i < m; i++) f[i] = f2cad::dint_(x+i); f2cad::Dependent(m, f); // declare f as dependent variable vector // check function value ok &= (f2cad::Value(f[0]) == 4.); ok &= (f2cad::Value(f[1]) == -3.); // Compute derivative of f w.r.t x using f2cad common interface. // This general interface works with adolc, cppad, and fadbad. // See GetStarted for package specific more efficient code. for(i = 0; i < m; i++) { double p = f2cad::Partial<doublereal>(i, i); ok &= f2cad::near_equal(p, 0., 1e-10, 1e-10); } if( ok ) return test_pass; return test_fail; } 
Input File: example/dint.cpp
6.12.1: dint.f
 double precision function dint (x) c***begin prologue dint c***revision october 1, 1980 c***category no. m2 c***keyword(s) truncation,greatest integer,double precision c***author fullerton w. (lasl) c***date written august 1979 c***purpose to find the largest integer whose magnitude does c not exceed x and convert to double precision. c***description c august 1979 edition. w. fullerton, los alamos scientific laboratory. c installed on the vax by dolores montano, c-3, 5/80. c c dint is the double precision equivalent of aint. this portable c version is quite efficient when the argument is reasonably small (a c common case), and so no faster machine-dependent version is needed. c c c***reference(s) c***routines called d1mach,i1mach,r1mach,xerror c***end prologue double precision x,xscl,scale,xbig,xmax,part,d1mach data npart,scale,xbig,xmax /0, 3*0.0d0 / c***first executable statement dint if(npart.ne.0) go to 10 ibase = i1mach(7) ndigd = i1mach(14) ndigi = min0(i1mach(8),i1mach(11)) - 1 nmax = 1.0d0/d1mach(4) xbig = amin1 (float(i1mach(9)),1.0/r1mach(4)) if (ibase.ne.i1mach(10)) call xerror ( 1 'dint algorithm error. integer base ne real base', 49, 2, 2) c npart = (ndigd + ndigi -1) / ndigi scale = ibase**ndigi c 10 if(x.lt.(-xbig) .or. x.gt.xbig) go to 20 c dint=int(sngl(x)) return c 20 xscl = dabs(x) if(xscl.gt.xmax) go to 50 c do 30 i=1,npart xscl=xscl/scale 30 continue c do 40 i=1,npart xscl = xscl*scale ipart=xscl part=ipart xscl=xscl - part dint = dint * scale + part 40 continue c if(x.lt.0.0d0) dint = -dint return c 50 call xerror ('dint dabs(x) may be too big to be represented 1as an exact integer', 67, 1, 1) dint=x return c end 
Input File: f.omh/dint.f.omh
6.13: dscal

6.13.a: Prototype
int f2cad::dscal_(integer *n, doublereal *da, doublereal *dx, integer *incx);

6.13.b: Fortran Source
6.13.1: dscal.f

6.13.c: Description
This example uses the routine dscal.f to compute   $\left( \begin{array}{c} f_0 \\ f_1 \end{array} \right) = a \left( \begin{array}{c} 1 \\ 2 \end{array} \right)$  Using the 8.3: f2cad_link routines f2cad::Independent and f2cad::Dependent, this defines the function   $f(a) = \left( \begin{array}{c} a \\ 2 a \end{array} \right)$  We check that the derivative of this function, calculated using the f2cad::Partial routine, satisfies   $f^{(1)} (a) = \left( \begin{array}{c} 1 \\ 2 \end{array} \right)$   # include <f2cad/dscal.hpp> test_result dscal(void) { bool ok = true; // Input values for dscal integer n = 1; doublereal a[1]; a[0] = 5.; // declare independent variables f2cad::Independent(n, a); integer m = 2; doublereal f[2]; f[0] = 1.; f[1] = 2.; integer incf = 1; // set f = a * f f2cad::dscal_(&n, a, f, &incf ); // declare dependent variables f2cad::Dependent(m, f); double p; p = f2cad::Partial<doublereal>(0, 0); ok &= f2cad::near_equal(p, 1., 1e-10, 1e-10); p = f2cad::Partial<doublereal>(1, 0); ok &= f2cad::near_equal(p, 2., 1e-10, 1e-10); if( ok ) return test_pass; return test_fail; } 
Input File: example/dscal.cpp
6.13.1: dscal.f
 subroutine dscal(n,da,dx,incx) c***begin prologue dscal c***date written 791001 (yymmdd) c***revision date 820801 (yymmdd) c***category no. d1a6 c***keywords blas,linear algebra,scale,vector c***author lawson, c. l., (jpl) c hanson, r. j., (snla) c kincaid, d. r., (u. of texas) c krogh, f. t., (jpl) c***purpose d.p. vector scale x = a*x c***description c c b l a s subprogram c description of parameters c c --input-- c n number of elements in input vector(s) c da double precision scale factor c dx double precision vector with n elements c incx storage spacing between elements of dx c c --output-- c dx double precision result (unchanged if n.le.0) c c replace double precision dx by double precision da*dx. c for i = 0 to n-1, replace dx(1+i*incx) with da * dx(1+i*incx) c***references lawson c.l., hanson r.j., kincaid d.r., krogh f.t., c *basic linear algebra subprograms for fortran usage*, c algorithm no. 539, transactions on mathematical c software, volume 5, number 3, september 1979, 308-323 c***routines called (none) c***end prologue dscal c double precision da,dx(1) c***first executable statement dscal if(n.le.0)return if(incx.eq.1)goto 20 c c code for increments not equal to 1. c ns = n*incx do 10 i = 1,ns,incx dx(i) = da*dx(i) 10 continue return c c code for increments equal to 1. c c c clean-up loop so remaining vector length is a multiple of 5. c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end 
Input File: omh/dscal.f.omh
6.14: idamax

6.14.a: Prototype
integer f2cad::idamax_(integer *n, doublereal *dx, integer *incx);

6.14.b: Fortran Source
6.14.1: idamax.f

6.14.c: Description
This example uses the routine idamax.f to determine the index corresponding to the maximum element in the vector   $\left( \begin{array}{cccc} 1 & 2 & 5 & 4 \end{array} \right)$  The maximum element is 5 and its index (in Fortran notation) is 3.  # include <f2cad/idamax.hpp> test_result idamax(void) { bool ok = true; // dx is a vector of length 4 double data[] = {1, 2, 5, 4}; doublereal dx[4]; integer i; for(i = 0; i < 4; i++) dx[i] = data[i]; // other arguments to idamax integer n = 4; integer incx = 1; // check return value ok &= ( 3 == f2cad::idamax_(&n, dx, &incx) ); if( ok ) return test_pass; return test_fail; } 
Input File: example/idamax.cpp
6.14.1: idamax.f
 integer function idamax(n,dx,incx) c***begin prologue idamax c***date written 791001 (yymmdd) c***revision date 820801 (yymmdd) c***category no. d1a2 c***keywords blas,double precision,linear algebra,maximum component, c vector c***author lawson, c. l., (jpl) c hanson, r. j., (snla) c kincaid, d. r., (u. of texas) c krogh, f. t., (jpl) c***purpose find largest component of d.p. vector c***description c c b l a s subprogram c description of parameters c c --input-- c n number of elements in input vector(s) c dx double precision vector with n elements c incx storage spacing between elements of dx c c --output-- c idamax smallest index (zero if n .le. 0) c c find smallest index of maximum magnitude of double precision dx. c idamax = first i, i = 1 to n, to minimize abs(dx(1-incx+i*incx) c***references lawson c.l., hanson r.j., kincaid d.r., krogh f.t., c *basic linear algebra subprograms for fortran usage*, c algorithm no. 539, transactions on mathematical c software, volume 5, number 3, september 1979, 308-323 c***routines called (none) c***end prologue idamax c double precision dx(1),dmax,xmag c***first executable statement idamax idamax = 0 if(n.le.0) return idamax = 1 if(n.le.1)return if(incx.eq.1)goto 20 c c code for increments not equal to 1. c dmax = dabs(dx(1)) ns = n*incx ii = 1 do 10 i = 1,ns,incx xmag = dabs(dx(i)) if(xmag.le.dmax) go to 5 idamax = ii dmax = xmag 5 ii = ii + 1 10 continue return c c code for increments equal to 1. c 20 dmax = dabs(dx(1)) do 30 i = 2,n xmag = dabs(dx(i)) if(xmag.le.dmax) go to 30 idamax = i dmax = xmag 30 continue return end 
Input File: omh/idamax.f.omh
6.15: initds

6.15.a: Prototype
integer f2cad::initds_(doublereal *dos, integer *nos, doublereal *eta);

6.15.b: Fortran Source
6.15.1: initds.f

6.15.c: Description
This example uses the routine initds.f to compute the minimum value of m such that   $\sum_{i=m}^{i=n} | a_i | \leq \eta$  Note that there are m terms in the rest of the summation; i.e.,   $\sum_{i=1}^{m-1} a_i$   # include <f2cad/initds.hpp> test_result initds(void) { bool ok = true; // a_0, a_1, a_2, a_3, a_4, a_5 doublereal a[] = { 1e0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5}; integer n = sizeof(a) / sizeof(a[0]); doublereal eta = 1e-3; integer m = f2cad::initds_(a, &n, &eta); ok &= (m == 4); if( ok ) return test_pass; return test_fail; } 
Input File: example/initds.cpp
6.15.1: initds.f
 function initds(dos,nos,eta) c***begin prologue initds c***date written 770601 (yymmdd) c***revision date 820801 (yymmdd) c***category no. c3a2 c***keywords chebyshev,double precision,initialize, c orthogonal polynomial,series,special function c***author fullerton, w., (lanl) c***purpose initializes the d.p. properly normalized orthogonal c polynomial series to determine the number of terms needed c for specific accuracy. c***description c c initialize the double precision orthogonal series dos so that initds c is the number of terms needed to insure the error is no larger than c eta. ordinarily eta will be chosen to be one-tenth machine precision c c input arguments -- c dos dble prec array of nos coefficients in an orthogonal series. c nos number of coefficients in dos. c eta requested accuracy of series. c***references (none) c***routines called xerror c***end prologue initds c double precision dos(nos) c***first executable statement initds if (nos.lt.1) call xerror ( 'initds number of coefficients lt 1', 1 35, 2, 2) c err = 0. do 10 ii=1,nos i = nos + 1 - ii err = err + abs(sngl(dos(i))) if (err.gt.eta) go to 20 10 continue c 20 if (i.eq.nos) call xerror ( 'initds eta may be too small', 28, 1 1, 2) initds = i c return end 
Input File: omh/initds.f.omh
6.16: xerror

6.16.a: Prototype
int f2cad::xerror_(char *messg, integer *nmessg, integer *nerr, integer *level, ftnlen messg_len);

6.16.b: Fortran Source
6.16.1: xerror.f

6.16.c: Description
This example uses the Cmlib routine xerror.f to print the message  Ok: xerror  (which is the expected output for this tests; see 7: run.cpp ).  # include <cstring> # include <f2cad/xerror.hpp> void xerror(void) { const char* messg = "Pass: xerror\n"; integer nmessg = std::strlen(messg); integer nerr = 1; integer level = 0; ftnlen messg_len = (ftnlen) nmessg; f2cad::xerror_(messg, &nmessg, &nerr, &level, messg_len); } 
Input File: example/xerror.cpp
6.16.1: xerror.f
 subroutine xerror(messg,nmessg,nerr,level) c***begin prologue xerror c***date written 790801 (yymmdd) c***revision date 820801 (yymmdd) c***category no. r3c c***keywords error,xerror package c***author jones, r. e., (snla) c***purpose processes an error (diagnostic) message. c***description c abstract c xerror processes a diagnostic message, in a manner c determined by the value of level and the current value c of the library error control flag, kontrl. c (see subroutine xsetf for details.) c c description of parameters c --input-- c messg - the hollerith message to be processed, containing c no more than 72 characters. c nmessg- the actual number of characters in messg. c nerr - the error number associated with this message. c nerr must not be zero. c level - error category. c =2 means this is an unconditionally fatal error. c =1 means this is a recoverable error. (i.e., it is c non-fatal if xsetf has been appropriately called.) c =0 means this is a warning message only. c =-1 means this is a warning message which is to be c printed at most once, regardless of how many c times this call is executed. c c examples c call xerror('smooth -- num was zero.',23,1,2) c call xerror('integ -- less than full accuracy achieved.', c 43,2,1) c call xerror('rooter -- actual zero of f found before interval f c 1ully collapsed.',65,3,0) c call xerror('exp -- underflows being set to zero.',39,1,-1) c c latest revision --- 19 mar 1980 c written by ron jones, with slatec common math library subcommittee c***references jones r.e., kahaner d.k., "xerror, the slatec error- c handling package", sand82-0800, sandia laboratories, c 1982. c***routines called xerrwv c***end prologue xerror character*(*) messg c***first executable statement xerror c call xerrwv(messg,nmessg,nerr,level,0,0,0,0,0.,0.) c statement above was commented out and the modified implementation c below is for standard linking to a c routine c c ------------------------------------------------------------------------- c This is a limited version of the specifications above for use with f2cad integer message(1000) n = min(nmessg, 1000) do i = 1, n message(i) = ichar( messg(i:i) ) 100 enddo call xerror2cpp(message, nmessg, nerr, level) return end 
Input File: omh/xerror.f.omh
7: Run Test of All the Library Routines

7.a: Syntax
example/adolc/run  example/cppad/run  example/fadbad/run

7.b: Description
These programs will test libf2c_adolc.a, libf2c_cppad.a, and libf2c_fadbad.a respectively. If all the test succeed, they will print  All the tests passed.  and return the value zero (no error). Otherwise they will print  At least one test failed.  and return a non-zero value (error).

7.c: Source Code
The source code for these programs follows below. Note that the definition of integer and doublereal, in the prototypes for the Fortran functions, depends on the preprocessor symbols (see 5: prototype ).  // system include files used for I/O # include <iostream> // define the types integer and doublereal # include <f2cad/f2cad.hpp> // ../add2lib.sh uses this comment for extern command extern test_result dbesks(void); extern test_result dbskes(void); extern test_result d9knus(void); extern test_result d9lgmc(void); extern test_result dgamma(void); extern test_result dint(void); extern test_result dcsevl(void); extern test_result dgamlm(void); extern test_result initds(void); extern test_result idamax(void); extern test_result daxpy(void); extern test_result ddot(void); extern test_result dgefa(void); extern test_result dgesl(void); extern test_result dscal(void); // xerror is a special case extern void xerror(void); // function that runs one test bool run(test_result test_fun(void), const char *name) { bool ok = true; switch( test_fun() ) { case test_none: std::cout << "None: " << name << std::endl; break; case test_pass: std::cout << "Pass: " << name << std::endl; break; case test_fail: std::cout << "Fail: " << name << std::endl; ok = false; break; } return ok; } // main program that runs all the examples int main(void) { bool ok = true; using std::cout; using std::endl; // ../add2lib.sh uses this comment for run command ok &= run( dbesks , "dbesks" ); ok &= run( dbskes , "dbskes" ); ok &= run( d9knus , "d9knus" ); ok &= run( d9lgmc , "d9lgmc" ); ok &= run( dgamma , "dgamma" ); ok &= run( dint , "dint" ); ok &= run( dcsevl , "dcsevl" ); ok &= run( dgamlm , "dgamlm" ); ok &= run( initds , "initds" ); ok &= run( idamax , "idamax" ); ok &= run( daxpy , "daxpy" ); ok &= run( ddot , "ddot" ); ok &= run( dgefa , "dgefa" ); ok &= run( dgesl , "dgesl" ); ok &= run( daxpy , "dscal" ); if( ok ) cout << endl << "Pass: None of the tests above failed." << endl; else cout << endl << "Fail: At least one test above failed." << endl; cout << " The next two lines should be identical:" << endl; // xerror is a special case because it is not checked automatically // and must be checked by eye. cout << "Pass: xerror" << endl; xerror(); cout << endl; if( ok ) return 0; return 1; } 
Input File: example/run.cpp
8: Utilities used to Test Library

8.a: Purpose
These library functions are used to test the f2cad library. They are not differentiable versions of original Fortran functions.

8.b: Contents

Input File: omh/utility.omh
8.1: Implementation of Cmlib i1mach Routine in Cpp

8.1.a: Syntax
j = i1mach_(i)

8.1.b: Purpose
This routine uses the C++ standard numeric_limits to implement the specifications required by the cmlib Fortran routine i1mach.

8.1.c: i
The argument i has prototype       integer* i  It specifies which machine constant to evaluate.

8.1.d: d
The return value i has prototype % integer %d% % and is the value of the requested machine constant.

8.1.e: Mapping
The meaning for the machine constants below was extracted from the comments in the Fortran routine i1mach. The constants corresponding to i = 1, 2, 3, 4 are not defined. If any of these constants are requested, i1mach will generate an assertion error. Also note that, in the f2cad library, single and double precision Fortran values are both represented using C++ doubles. Hence, the return values for single and double precision are the same.
 i j 1 standard input unit. 2 standard output unit. 3 standard punch unit. 4 standard error message unit. 5 number of bits per integer storage unit. 6 number of characters per integer storage unit. 7 base used for representation of integers. 8 number of base digits per integer storage unit. 9 largest integer magnitude 10 base for representing floating point numbers. 11 number of base digits in representation of single precision. 12 smallest exponent in representation of single precision. 13 largest exponent in representation of single precision. 14 number of base digits in representation of double precision. 15 smallest exponent in representation of double precision. 16 largest exponent in representation of double precision.

Input File: lib/i1mach.cpp
8.2: Implementation of Cmlib d1mach Routine in Cpp

8.2.a: Syntax
d = d1mach_(i)

8.2.b: Purpose
This routine uses the C++ standard numeric_limits to implement the specifications required by the cmlib Fortran routine d1mach.

8.2.c: i
The argument i has prototype       integer* i  It specifies which machine constant to evaluate.

8.2.d: d
The return value d has prototype % doublereal %d% % and is the value of the requested machine constant.

8.2.e: Mapping
The Fortran definition for the machine constants below was extracted from the comments in the Fortran routine d1mach. The C++ definitions were extracted from the numeric_limits section of the C++ standard.
 i d 1 B**(EMIN-1), the smallest positive magnitude. 2 B**EMAX*(1 - B**(-T)), the largest magnitude. 3 B**(-T), the smallest relative spacing. 4 B**(1-T), the largest relative spacing. 5 LOG10(B).

Input File: lib/d1mach.cpp

8.3.a: Syntax
void f2cad::Independent(n, x)  void f2cad::Dependent(m, y)  double f2cad::Partial<doublereal>(i, j)  double f2cad::Value(z)  double f2cad::Sign(x, y) 
8.3.b: Purpose

8.3.c: doublereal
The type doublereal is the C++ type that corresponds to Fortran single and double precision. This type is mapped to an AD type that depends on the library being used:
 Library typedef libf2c_adolc.a typedef adouble         doublereal libf2c_cppad.a typedef CppAD:: doublereal libf2c_fadbad.a typedef F       doublereal

8.3.d: integer
The type integer is the C++ type that corresponds to Fortran integers.

8.3.e: Vector Elements

8.3.f: Tape State
Each doublereal type has a tape state and the initial tape state is Empty.

8.3.g: Independent
The argument n has prototype       integer n  It specifies the dimension of the domain space for the function being defined. The argument x has prototype       doublereal *x  It is a vector of length n that contains the independent variables (and values) for this function recording. The value of the elements of x are not affected by calling this routine. (The elements are not const in the prototype because the AD information connected to them may be affected.)

8.3.g.a: Tape State
The tape state must be either Empty or Complete when this routine is called. After this routine is called, the tape state will be Recording,

8.3.h: Dependent
The argument m has prototype       integer m  It specifies the dimension of the range space for the function being defined. The argument y has prototype       doublereal *x  It is a vector of length m that contains the dependent variables (and values) for this function recording. The value of the elements of y are not affected by calling this routine. (The elements are not const in the prototype because the AD information connected to them may be affected.)

8.3.h.a: Tape State
The tape state must be Recording when this routine is called. After this routine is called, the tape state will be Complete,

8.3.i: Partial
The argument i has prototype       integer i  It must be greater than or equal zero and less than the dimension of the range space for the function. The argument j has prototype       integer j  It must be greater than or equal zero and less than the dimension of the domain space for the function. The return value from Partial is the partial derivative of the i-th dependent variable with respect to the j-th dependent variable (i-th function value with respect to j-th function argument).

8.3.i.a: Tape State
The tape state must be Complete when this routine is called. The tape state is not affected by this call.

8.3.j: Value
The argument z has prototype       const doublereal &z  This function returns the value corresponding to the doublereal object z .

8.3.k: Sign
The arguments x and y have prototypes       doublereal* x      doublereal* y  This function returns   $\left\{ \begin{array}{ll} + |x| & \R{if} \; y \geq 0 \\ - |x| & \R{otherwise} \end{array} \right.$ 

8.3.k.a: Tape State
The tape state must be either Empty or Complete when this routine is called. The tape state is not affected by this call.

8.3.l: Source Code
The following table contains links to the source code for each AD type:

8.3.m: Contents

 // assume that F2CAD_USE_Adolc is defined # include <f2cad/f2cad.hpp> namespace { static bool Recording = false; // currently recording a function static int Domain = 0; // domain dimension static int Range = 0; // range dimension static double *Xmemory = 0; // domain value and differential static double **X = 0; static double *Ymemory = 0; // range value and differential static double **Y = 0; } namespace f2cad { void Independent(integer n, doublereal *x) { // check not two Independents in a row assert( ! Recording ); // make sure the set of independent variables not empty assert( n > 0 ); // check for previously defined function (non-zero pointers) if( Xmemory != 0 ) { assert( Xmemory != 0 ); delete [] Xmemory; assert( Ymemory != 0 ); delete [] Ymemory; assert( X != 0 ); delete [] X; assert( Y != 0 ); delete [] Y; } // allocate memory for independent variable value and differential // and store domain value in first column of X Xmemory = new double[2 * n]; X = new double *[n]; int j; for(j = 0; j < n; j++) { X[j] = Xmemory + 2 * j; X[j][0] = x[j].value(); } // start recording on tape number 1 Recording = true; Domain = n; Range = 0; int tag = 1; int keep = 1; trace_on(tag, keep); // declare independent variables for(j = 0; j < n; j++) x[j] <<= X[j][0]; return; } void Dependent(integer m, doublereal *y) { // check that directly follows call to Independent assert( Recording ); assert( Range == 0 ); // make sure the set of dependent variables is not empty assert( m > 0 ); // declare the independent variables int i; double yi; for(i = 0; i < m; i++) y[i] >>= yi; // allocate memory for dependent variable value and differential // and store domain value in first column of X Ymemory = new double[2 * m]; Y = new double *[m]; for(i = 0; i < m; i++) Y[i] = Ymemory + 2 * i; // stop the recording Range = m; Recording = false; trace_off(); } // linked here from Partial<doublereal>(integer i, integer j) double Partial(const doublereal &zero, integer i, integer j) { // check that function is defined assert( ! Recording ); assert( Domain > 0 ); // check dimension of arguments int m = Range; int n = Domain; assert( 0 <= i ); assert( 0 <= j ); assert( i < m ); assert( j < n ); assert( zero == doublereal(0) ); // move data into structure expected by Adolc int k; for(k = 0; k < Domain; k++) X[k][1] = 0.; X[j][1] = 1.; // compute dy int tag = 1; // tape identifier int keep = 1; // keep data in tape int d = 1; // highest derivative degree forward(tag, m, n, d, keep, X, Y); // return the result return Y[i][1]; } double Value(const doublereal &z) { doublereal tmp = z; return tmp.value(); } doublereal Sign(const doublereal* x, const doublereal* y) { if( *y >= 0. ) return abs(*x); return -abs(*x); } } // end f2cad namespace 
 // assume that F2CAD_USE_CppAD is defined # include <f2cad/f2cad.hpp> namespace { static CppADvector< CppAD::AD<double> > X; static CppAD::ADFun<double> *F = 0; } namespace f2cad { void Independent(integer n, doublereal *x) { // check not two Independents in a row assert( X.size() == 0 ); // make sure the set of independent variables not empty assert( n > 0 ); // check for a previously defined function (nonzero pointer) if( F != 0 ) { delete F; F = 0; } // move data into structure expected by CppAD X.resize(n); integer j; for(j = 0; j < n; j++) X[j] = x[j]; // declare independent variables CppAD::Independent(X); // make x equivalent to X for(j = 0; j < n; j++) x[j] = X[j]; return; } void Dependent(integer m, doublereal *y) { // check that directly follows call to Independent assert( X.size() > 0 ); assert( F == 0 ); // make sure the set of dependent variables not empty assert( m > 0 ); // move data into structure expected by CppAD CppADvector< CppAD::AD<double> > Y(m); integer i; for(i = 0; i < m; i++) Y[i] = y[i]; // construct the AD function object F = new CppAD::ADFun<double>(X, Y); // no longer need independent variable information X.resize(0); } // linked here from Partial<doublereal>(integer i, integer j) double Partial(const doublereal &zero, integer i, integer j) { // check that directly follows call to Dependent assert( X.size() == 0 ); assert( F != 0 ); // check dimension of arguments integer m = F->Range(); integer n = F->Domain(); assert( 0 <= i ); assert( 0 <= j ); assert( i < m ); assert( j < n ); assert( zero == doublereal(0) ); // move data into structure expected by CppAD CppADvector<double> dx(n); integer k; for(k = 0; k < n; k++) dx[k] = 0.; dx[j] = 1.; // Forward mode calculation CppADvector<double> dy(m); dy = F->Forward(1, dx); // return the result return dy[i]; } double Value(const doublereal &z) { // convert in a way that works even if we are recording return CppAD::Value( Var2Par(z) ); } doublereal Sign(const doublereal* x, const doublereal* y) { if( *y >= 0. ) return abs(*x); return -abs(*x); } } // end f2cad namespace 
 // assume that F2CAD_USE_Fadbad is defined # include <f2cad/f2cad.hpp> namespace { static bool Recording = false; // currently recording a function static int Domain = 0; // domain dimension static int Range = 0; // range dimension static doublereal *Fun = 0; // currently defined function object } namespace f2cad { void Independent(integer n, doublereal *x) { // check not two Independents in a row assert( ! Recording ); // make sure the set of independent variables not empty assert( n > 0 ); // delete old function object if( Fun != 0 ) { delete [] Fun; Fun = 0; } // declare independent variables int j; for(j = 0; j < n; j++) x[j].diff(j, n); // change to recording state Recording = true; Domain = n; Range = 0; return; } void Dependent(integer m, doublereal *y) { // check that directly follows call to Independent assert( Recording ); assert( Range == 0 ); assert( Fun == 0 ); // make sure the set of dependent variables is not empty assert( m > 0 ); // create the function object Fun = new doublereal[m]; int i; for(i = 0; i < m; i++) Fun[i] = y[i]; // stop the recording Range = m; Recording = false; } // linked here from Partial<doublereal>(integer i, integer j) double Partial(const doublereal &zero, integer i, integer j) { // check that function is defined assert( ! Recording ); assert( Fun != 0 ); assert( 0 <= i ); assert( 0 <= j ); assert( i < Range ); assert( j < Domain ); assert( zero == doublereal(0) ); // compute dy return Fun[i].d(j); } double Value(const doublereal &z) { // value member function is not const doublereal tmp = z; return tmp.x(); } doublereal Sign(const doublereal* x, const doublereal* y) { if( *y >= 0. ) return abs(*x); return -abs(*x); } } // end f2cad namespace 
8.4: Determine if Two Double Values Are Nearly Equal
 Syntax # include  bool near_equal(double x, double y, double r, double a)

8.4.a: Description
Returns true if   $x = y \; {\rm or} \; | x - y | \leq a \; {\rm or} \; \frac{ | x - y | } { |x| + |y| } \leq r$  Otherwise, it return false. Note that r is a relative error bound and a is an absolute error bound.

Each of the libraries libf2c_adolc.a, libf2c_cppad.a, and libf2c_fadbad.a has a copy of this routine in it.
Input File: lib/near_equal.cpp
8.5: Implementation of Cmlib xerror Routine in Cpp

8.5.a: Syntax
xerror2cpp_(message, n, nerr, level)

8.5.b: Description
The f2cad version of 6.16.1: xerror.f uses this routine for its implementation. This is not completely faithful to the specifications for xerror.

8.5.c: message
The argument message has prototype       const char *message  It points to the beginning of the error message to be printed.

8.5.d: n
The argument n has prototype       const int *n  It is the number of characters in the error message.

8.5.e: nerr
The argument nerr has prototype       const int *nerr  It is a number that identifies the error.

8.5.f: level
The argument level has prototype       const int *level  It has the following meaning:
 -1 a warning, message is printed (but not multiple times in a row) 0 a warning, message is printed 1 a non-fatal error, error number and message are printed 2 a fatal error, error number and message are printed and program aborts

Input File: lib/xerror2cpp.cpp

9.a: Required Software

9.b: Example Steps
1. Convert to lower case and copy the Fortran source code for the routine you are adding to the fortran subdirectory of the distribution directory. For the mysub example, execute the following command in the distribution directory:       tr 'A-Z' 'a-z' < mysub.f > fortran/mysub.f 
2. In the f2cad distribution directory, execute the command       ./add2lib.sh name  where name is the name of the Fortran routine. For the mysub example, execute the following command       ./add2lib.sh mysub 
3. In the f2cad distribution directory, execute the command       ./add2lib.sh use_new_files  This command will copy the *.new files to the corresponding files (but first it will back up the corresponding files in *.old ). If name is not mysub, the initial version of       example/name.cpp  will return test_none. You can edit this file to create a real test and example of name .
4. In the file build.sh, make sure that F2C_DIR, ADOLC_DIR, and FADBAD_DIR have the proper values for your system. Then execute the command  ./build.sh all 
5. You can undo the effect of the add2lib.sh commands above by executing  ./add2lib.sh use_old_files 

10.a: Introduction
This section contains a list of changes to f2cad in reverse order by date. Its purpose is to assist you in learning about changes between various versions.

10.b: 10-04-28
1. Added the routines 6.5: dbskes and 6.4: dbesks .
2. Bring install instructions up to date.

10.c: 10-04-27
1. Added the routine 6.1: d9knus .
2. Added instructions for converting Fortran to lower case to 9: add2lib.sh .
3. Improved sed script that converts from C to C++ so that it handles more cases.
4. Fix undefined external MAIN_ on Fedora where using yum installed version of f2c (related to dynamic link -lf2c).

10.d: 10-04-26
1. Update 3: install for move to CoinBazaar.
2. Change the configure flag CPP_ERROR_WARN to CXX_FLAG.

10.e: 10-04-25
1. Change the test return type to  typedef enum { test_none, test_pass, test_fail } test_return;  This enables the driver program to inform the user when there is no test for a routine (which is how 9: add2lib.sh initializes each example / test).
2. Added the 8.1: i1mach utility.
3. Added the 6.12: dint , 6.2: d9lgmc , and 6.9: dgamma library functions.

10.f: 10-04-24
Alphabetize the 6: library section. Put a direct link to the Fortran source code directly after the prototype for each subsection of 6: library (the Fortran source has the documentation).

10.g: 10-04-04
1. Move f2cad source code to https://projects.coin-or.org/CoinBazaar/browser
2. Fix an upper/lower case bug that did not show up under cygwin (which is case insensitive on file names).
3. Remove some unused files.

10.h: 10-04-03
1. Moved example/*.h to example/*.cpp .
2. Separated the example build for each packages into its own directory. For example, moved example/RunAdolc to example/adolc/adolc, example/RunCppAD to example/cppad/cppad, and example/RunFadbad to example/fadbad/fadbad.
3. Change C++ include file names from *.h to *.hpp .
4. Fix mistake in function all return value prototypes (stored in f2cad/*.omh ).
5. Make fortran/mysub.f a working example for using 9: add2lib.sh .

10.i: 10-04-02
Added the routine 6.6: dcsevl .

10.j: 10-03-31
1. Convert Blas routine names to the standard for double precision; e.g., saxpy was changes to 6.3: daxpy .
2. Added the routines 8.2: d1mach and 6.8: dgamlm .

10.k: 10-03-29
1. Modifications to bring up to date so worked with current versions of all the tools used.
2. Fix a bug in 4.1: start_adolc.cpp (accessing memory after it had been freed).
3. Include original Fortran and add links in examples to: 6.14.1: idamax.f , 6.3.1: daxpy.f , 6.7.1: ddot.f , 6.10.1: dgefa.f , 6.11.1: dgesl.f , 6.13.1: dscal.f .
4. Add the routines 6.16: xerror and 8.5: xerror2cpp to the library.
5. Add the cmlib routine 6.15: initds .

10.l: 09-03-01
The location of f2c on the web changed and its instructions no longer worked. An f2c_install.sh script was written which automates most of the install process (and works with the new location of f2c on the web).  Upgrade to new versions of 3.h: Adolc , and 3.j: Fadbad .

10.m: 05-08-27
The 9: add2lib.sh example fortran/mysub.f has been added. This shows to add multiple fortran functions and subroutines in one file. The advantage (or disadvantage) of this is that you only have to test one of the routines in the file.

10.n: 05-08-25
A script file that makes it easier to add new fortran routines to the libraries 9: add2lib.sh was added.

10.o: 05-04-09
The first version of f2cad was uploaded to        http://www.seanet.com/~bradbell/f2cad.unix.tar.gz .

10.p: 05-04-10
Use the vector   f` for the function value in the example and test of 6.3: daxpy , 6.7: ddot , and 6.13: dscal .
Input File: omh/whats_new.omh
11: Alphabetic Listing of Cross Reference Tags

12: Keyword Index
A
absolute
difference 8.4: Determine if Two Double Values Are Nearly Equal
all
run examples 7: Run Test of All the Library Routines
B
blas
daxpy 6.3: daxpy
ddot 6.7: ddot
dgefa 6.10: dgefa
dgesl 6.11: dgesl
dscal 6.13: dscal
idamax 6.14: idamax
C
C++
d1mach 8.2: Implementation of Cmlib d1mach Routine in Cpp
i1mach 8.1: Implementation of Cmlib i1mach Routine in Cpp
cmlib
initds 6.15: initds
xerror 6.16: xerror
D
d1mach
C++ 8.2: Implementation of Cmlib d1mach Routine in Cpp
d9knus 6.1: d9knus
d9lgmc 6.2: d9lgmc
daxpy
blas 6.3: daxpy
dbesks 6.4: dbesks
dbskes 6.5: dbskes
dcsevl
fnlib 6.6: dcsevl
ddot
blas 6.7: ddot
derivative
dgamlm
fnlib 6.8: dgamlm
dgamma 6.9: dgamma
dgefa
blas 6.10: dgefa
dgesl 6.10: dgefa
blas 6.11: dgesl
difference
absolute or relative 8.4: Determine if Two Double Values Are Nearly Equal
dint 6.12: dint
double
near equal 8.4: Determine if Two Double Values Are Nearly Equal
dscal
blas 6.13: dscal
E
equal
near 8.4: Determine if Two Double Values Are Nearly Equal
error
xerror 8.5: Implementation of Cmlib xerror Routine in Cpp
example
run 7: Run Test of All the Library Routines
F
fnlib
dcsevl 6.6: dcsevl
dgamlm 6.8: dgamlm
fortran
prototype 5: C++ Types Corresponding to Fortran Prototypes
G
getting
I
i1mach
C++ 8.1: Implementation of Cmlib i1mach Routine in Cpp
idamax
blas 6.14: idamax
initds
cmlib 6.15: initds
interface
L
P
prototype
fortran 5: C++ Types Corresponding to Fortran Prototypes
R
relative
difference 8.4: Determine if Two Double Values Are Nearly Equal
run
examples 7: Run Test of All the Library Routines
S
start