Showing posts with label Fortran. Show all posts
Showing posts with label Fortran. Show all posts

Wednesday, November 30, 2011

NOTES GSL

Prerequisites
Building using Visual Studio C/C++
GSL Examples
GSL Considerations
Build GSL with MinGW
Running MSVC with minGW-built-GSL
GSL Makefile
General C define macros and typedef
GSL-C-Fortran Framework

Prerequisites
===============
For running on Windows system:

1. mingw (minimalist GNU for Windows)
- go to http://www.mingw.org/
- Download and install both mgw and msys: eg. mingw-get-inst-20110802.exe
- Install and select C++, Fortran and MSYS options

2. gsl (Gnu Scientific Library)
- go to http://www.gnu.org/s/gsl/
- Download from nearest GNU mirror, eg. http://ftpmirror.gnu.org/gsl/.
- Unpack it to a general place like c:\gsl....

Building using Visual Studio C/C++
===================================
http://gladman.plushost.co.uk/oldsite/computing/gnu_scientific_library.php
http://www.quantcode.com/modules/smartfaq/faq.php?faqid=33


The GSL subroutines are grouped according to functionality. For this document, we say that these functionalities are grouped into packages. Each package is a folder containing the source code for that functionality. In this section, we will compile the package called BLAS.

Note that this BLAS is actually gls_blas. It uses the underlying blas implementation called cblas.
cblas is the standard blas for C/C++ and is also available in MKL.
In this example build, we will build the gsl_blas and configure is to use cblas from MKL.

1. Create a Visual Studio C/C++ Project called GSLProj, under the Solution called GSLSoln.
Add New Project - Visual C++ - Win32 - Win32 Project - {DLL, Empty Project}

2. Go to Windows Explorer and copy the blas folder from c:\gsl....
and paste into the GSLProj folder.
Update: The following packages seemed to be a bare minimum and have enabled GSL to be compiled successfully with MSVC 2008 - blas, block, err, ieee-utils, matrix, sys, test, vector.
Update: Second batch of packages to be added:
- cdf, cheb, complex, eigen, integration, linalg, min, multimin, multiroots, permutation, poly, randist, rng, roots, sort, specfun,

3. In Windows Explorer, copy these files from c:\gsl
    config.h , templates_on.h, templates_off.h, build.h
and paste into the GSLProj folder.

4. In Windows Explorer, under c:\gsl...., do a search for gsl_*.h. This will find all the gsl header files in the source code. Copy all of these and paste into a folder called gsl and put this folder under the GSLProj folder.

5. Go back to the Visual Studio GSLProj, go to the Solution Explorer and click on the little icon called "Show All Files" just on top of the Solution tree.

6. In the Solution Explorer, the tree becomes a folder tree view. Go to the branch GSLProj, right click on the blas folder (which appears now) and select "Include in Project". Inlcude any other project as necessary, only include their *.c is OK. All header files do not need to be INCLUDED via Visual Studio because the C compiler will pull the contents of *.h hearder files in.

7. For the ieee-utils package, exclude the files below which are for various different operating systems other than Windows. These files have names like: fp-[OS].c
Including fp.c seems to be OK.

For the matrix package, Include the matrix folder into the project.
Then exclude files: *.h, *.lo, *_source.c
Ensure that normal c files are included:   *.c
Note that the *_source.c files are actually pulled into the base *.c file similar to how header files are included.

Specfun package also has these *.c files to be excluded from compile:
cheb_eval.c, cheb_eval_mode.c
 
In many matrix files, eg copy_source.c, file_source.c, getset_source.c, the copy.c, file.c, getset.c uses the templates_on/off.h which dynamically creates code using Preprocessor Macro definitions. This causes SEVERE problems at the line:
#define SHORT complex
because Microsoft Visual C redefines this in their math.h file as:
#define complex _complex
DO  NOT Put this macro definition /D__cplusplus  which may solve the _complex problem but creates new problems with other MS files like:

Error 1 error C2061: syntax error : identifier 'vc_attributes' c:\program files\microsoft visual studio 9.0\vc\include\codeanalysis\sourceannotations.h 42 VCPP_testDLL01

To solve the _complex problem definition, add this at the top of templates_on.h:
#ifdef complex
#define fixCOMPLEX_ 1
#undef complex
#endif

Also add this at the end of templates_off.h
// ASSUMES that the predefined complex is _complex in Microsoft's VC include\math.h and crt\src\math.h
#if fixCOMPLEX_ == 1
#define complex _complex
#undef fixCOMPLEX_
#endif

8. Right click on the GSLProj icon and select properties. In the Properties configuration, Configuration Properties:
- under C/C++ - Command Line
Add these compiler options: /DGSL_DLL /DDLL_EXPORT
- under C/C++ - General - Additional Include Directories:
      enter the path to the current project folder.
 enter the path to \mkl\include
- under Linker - General - Additional Library Directories:
      enter \fortran\lib\ia32
 enter \mkl\ia32\lib
- under Linker - Input - Additional Dependencies, enter these mkl, openmp libraries
   for 32bit:
mkl_intel_c_dll.lib
mkl_intel_thread_dll.lib
mkl_core_dll.lib
libiomp5md.lib
   for 64bit:
mkl_intel_lp64_dll.lib
mkl_intel_thread_dll.lib
mkl_core_dll.lib
libiomp5md.lib

The setup above enable the project to be linked to MKL's cblas routines.

9. Edit the config file so that the "!" is removed in the following lines:
#if HAVE_DECL_LDEXP
#if HAVE_DECL_FREXP
#if HAVE_DECL_HYPOT
For some reason the use of preprocessor macros does not seem to work well. Hence the config.h file need to be modified as above.

10. Removing cblas from GSL code.
- Exclude the cblas directory from the project
- Exclude the file gsl\gsl_cblas.h
- Remove the line
     #include
  and replace if necessary with
     #include
  for these files:
blas\blas.c
doc\examples\cblas.c    // probably not needed
eigen\francis.c
eigen\nonsymmv.c
gsl_blas_types.h

11. Due to using mkl_cblas.h, these definitions need to be changed.
in gsl_blas_types.h, change to:
typedef  CBLAS_ORDER       CBLAS_ORDER_t;
typedef  CBLAS_TRANSPOSE   CBLAS_TRANSPOSE_t;
typedef  CBLAS_UPLO        CBLAS_UPLO_t;
typedef  CBLAS_DIAG        CBLAS_DIAG_t;
typedef  CBLAS_SIDE        CBLAS_SIDE_t;
... the enum keyword has been removed.

12. Add the following Preprocessor definitions
/P    - only to check the preprocessed intermediate files

13. To add cdf package:
- add this to beta_inc.c
#include
- exclude from project, usually because these files are included in other files and are not really standalone files:
test_auto.c, beta_inc.c

cdf package needs the following packages:
specfun, cheb


14. Various packages of GSL, such as block and matrix, have some files with the same names such as file.c.
But the way MSVC compiles and builds by putting all *.obj files into the same DEBUG or RELEASE folder causes a problem. A simple way is to rename some of the clashing file names. This does not affect the ability to use the GSL-DLL functions, because functions have different names withing files with same names such as file.c

To quickly identify where duplicate filenames exist, go to Solution Explorer and click on the GSLProj project folder.
Near the top of the Solution Explorer panel, click on the icon to switch to Show All Files, such that all *.c files appear under the Source Files folder under your GSLProj project folder. Then rename any duplicate file names *.c.

Alternative - Right-click, Properties, C++, Output Files. Change the Object File Name from $(IntDir)\ to, say, $(IntDir)\block\. You can repeat for the vector file group, etc.

15. Inline functions.
The following packages are affected: cdf
MSVC++ accepts both inline and _inline BUT when the code is *.c and compiled with MSVC, then it only recognize _inline.
Error Messages are like:
Error 1 error C2054: expected '(' to follow 'inline' ...\specfunc\cheb_eval_mode.c
Error 2 error C2085: 'cheb_eval_mode_e' : not in formal parameter list ...\specfunc\cheb_eval_mode.c
Error 3 error C2143: syntax error : missing ';' before '{' ...\specfunc\cheb_eval_mode.c 6
Solution:
In specfun/cheb_eval_mode.c, airy.c , gamma.c, bessel.c, bessel_olver, bessel_zero, coupling.c, ellint.c
erfc.c, legendre_con, trig.c, zeta.c
change from inline to _inline

16. Notes on cdf and specfun packages.
- see 15. for "inline" issue
- Files to be renamed from *.c to *_cdf.c are:
     cdf/beta.c   to avoid clash with specfun/beta.c
- files to be excluded from compile:
cdf/beta_inc.c

 17. Add randist package - needed by cdf
 - inline keyword modified, see note 15; for files:
     binomial_tpe.c, discrete.c, shuffle.c
 - Files to be renamed from *.c to *_randist.c are:
     beta.c, binomial.c, cauchy.c, chisq.c, exponential.c, exppow.c, fdist.c, flat.c, gamma.c,
gauss.c, geometric.c, hyperg.c, laplace.c, logistic.c, lognormal.c, nbinomial.c, pareto.c, pascal.c, poisson.c, rayleigh.c, rdist.c, weibull.c

18. Add rng package - needed by randist
 - inline keyword modified, see note 15; for files:
     ... almost all files ...
 - Files to be renamed from *.c to *_rng.c are:
     file.c

19. Add complex package - needed by specfun
 - Files to be renamed from *.c to *_complex.c are:        
     inline.c

20. Add eigen package - needed by specfun
 - inline keyword modified, see note 15; for files:
     gen.c, francis.c, qrstep.c, jacobi.c,
- Remove           //#include
  and replace by     #include
  in these files:    francis.c, nonsymmv.c
- exclude from project:
      qrstep.c


21. Add linalg package - needed by
 - inline keyword modified, see note 15; for files:
     apply_givens.c, cholesky.c, givens.c, exponential.c
- exclude from project:
     apply_givens.c, givens.c, svdstep.c


22. Add permutation package - needed by linalg
 - inline keyword modified, see note 15; for files:
 - Files to be renamed from *.c to *_permutation.c are:        
      file.c, init.c, inline.c,
- exclude from project:
      permute_source.c


Add integration package - needed by randist
 - inline keyword modified, see note 15; for files:
    append.c, initialise.c, positivity.c, set_initial.c, qpsrt.c, util.c, reset.c, qelg.c,
 - Files to be renamed from *.c to *_integration .c are:        
     
- exclude from project:
     append.c, cquad_const.c, err.c, initialise.c, positivity.c, ptsort.c, qpsrt.c, util.c, reset.c
qelg.c, qc25c.c, qc25f.c, qc25s.c, qpsrt2.c, set_initial.c,

NOTE: if some *.c files which are included by other *.c files are NOT excluded from compile, then errors like these occur:
Error 1 error C2143: syntax error : missing ')' before '*'
Error 2 error C2143: syntax error : missing '{' before '*'
Error 3 error C2059: syntax error : 'type'
Error 4 error C2059: syntax error : ')'


Add sort package - needed by eigen
 - inline keyword modified, see note 15; for files:
     sort.c, sortind.c, sortvec_source.c, sortvecind_source.c,
 - Files to be renamed from *.c to *_sort.c are:        
     
- exclude from project:
     sortvec_source.c, sortvecind_source.c, subset_source.c,  subsetind_source.c,
test_heapsort.c, test_source.c,


Add min package - for minimization routines
-  Add static to these functions:
      min/test.c::my_error_handler()


Add roots package
 - inline keyword modified, see note 15; for files:
 - Files to be renamed from *.c to *_roots.c are:        
      test.c, test_funcs.c
- exclude from project:


Add multimin package
 - inline keyword modified, see note 15; for files:
      simplex2.c
 - Files to be renamed from *.c to *_multimin.c are:
      test.c, test_funcs.c
- exclude from project:
directional_minimize.c, linear_minimize.c, linear_wrapper.c,


Add poly package - needed by multimin
 - inline keyword modified, see note 15; for files:
 - Files to be renamed from *.c to *_multimin.c are:
      test.c, eval.c, deriv.c,
- exclude from project:
 companion.c, balance.c, qr.c,


Add multiroots package
 - inline keyword modified, see note 15; for files:
 - Files to be renamed from *.c to *_multiroots.c are:
      convergence.c, fdfsolver.c, fsolver, newton.c, test_funcs.c, test.c
- exclude from project:
      enorm.c, dogleg.c




23. To make use of GSL's own test, do the following to the selected GSL packages as needed:
cdf, specfunc, rng, randist, min, roots, multimin, poly,

- eg. in cdf/test.c  (which has been renamed to test_cdf.c), rename the function
     main      to    main_cdf.c

- ensure say, cdf/test.c  is now included in the project and compile.

- for the test functions, say: void main_cdf(void)
  put its prototype function in your own driver program using extern, like:
     extern "C" { void main_cdf(void);  }

- in the main functions, say main_cdf(), at the last line of the function, change this:
     exit(gsl_test_summary());
  to this, to return the EXIT_SUCCESS or EXIT_FAILURE code
     return gsl_test_summary();

- in the main functions, say main_sf(), if there are arguments that are not used like this:
     int main_sf(int argc, char * argv[])
  to the following
     int main_sf(void)

- in cdf package ONLY: in file test_cdf.c, change the line from:
     void test_gamma (void)
  to rename it to:
     void test_gamma_cdf (void)
  and in the same file, change the call to the test_gamma function from:
     test_gamma();
  to the new name :
     test_gamma_cdf();
  This is done to prevent clash with a function of the same name test_gamma() in specfunc.

- in randist package only, change the following functions by adding the keyword "static",
  so that other files to not uses these functions. The functions made "static" are:
     test_beta, test_binomial_pdf, test_chisq, test_exponential, test_fdist, test_gamma,
test_ugaussian, test_geometric_pdf, test_hypergeometric2_pdf, test_negative_binomial_pdf,
test_pascal_pdf, test_poisson_pdf

- in roots package only, change the following functions by adding the keyword "static",
  so that other files to not uses these functions. The functions / variables made "static" are:
     test_roots.c: EPSREL, EPSABS, MAX_ITERATIONS, test_f, test_f_e, my_error_handler,
test_funcs_roots.c: create_function, func1, func2, func3, func4

- in multimin package only, change the following functions by adding the keyword "static",
  so that other files to not uses these functions. The functions / variables made "static" are:
     test_fdf, test_f,

- in multiroots package only, change the following functions by adding the keyword "static",
  so that other files to not uses these functions. The functions / variables made "static" are:
     test_multiroots.c: test_fdf, test_f,
     test_funcs_multiroots.c: rosenbrock, rosenbrock_initpt, rosenbrock_f, rosenbrock_df, rosenbrock_fdf,
wood, wood_initpt, wood_f, wood_df, wood_fdf,
     roth, roth_initpt, roth_f, roth_df, roth_fdf,

- Strange reason, the two functions gsl_cdf_logistic_Q and  gsl_cdf_logistic_P are not visible to the test program and cannot be compiled. Typical error messages are:
Error 1  error LNK2001: unresolved external symbol gsl_cdf_logistic_Q test_cdf.obj vcpp_gsl
Error 1  error LNK2019: unresolved external symbol gsl_cdf_logistic_Q test_cdf.obj vcpp_gsl
    to solve this problem, go to the file cdf/logistic.c and copy the two functions of
gsl_cdf_logistic_Q and  gsl_cdf_logistic_P
    and paste it to cdf/test_auto.c and place it above the line:
void test_auto_logistic (void)
    ALTERNATIVE: Better solution is to check no two source files has the same name. This is already described in point 15. When two files have same name, one of them is not compiled and so its functions cannot be found.

- other files where functions are already identified somewhere else in test*.c. The solution is to add the static keyword to these functions:
   test_integration.c: my_error_handler
   test_linalg.c: my_error_handler


Debugging:
Modified test_roots.c: void test_f, test\results.c: void gsl_test()
Completed [12017/12017]
GSL cdf test result is 0
 end of gsl cdf test
Completed [14995/14995]
GSL specfunc test result is 0
 end of gsl specfunc test
-10, 1072693247, 0, 0
FAIL: results.c::gsl_test:: [15003]incorrect precision (%g obs vs %g expected)
-16, 1072693247, -350469331, 1058682594
FAIL: results.c::gsl_test:: [15005]incorrect precision (%g obs vs %g expected)
2142418408, 1061077358, 0, 0
FAIL: results.c::gsl_test:: [15007]incorrect precision (%g obs vs %g expected)
-13, 1072693247, -350469331, 1058682594
FAIL: results.c::gsl_test:: [15020]incorrect precision (%g obs vs %g expected)
2142401953, 1061077358, 0, 0
FAIL: results.c::gsl_test:: [15022]incorrect precision (%g obs vs %g expected)
-10, 1072693247, 0, 0
FAIL: results.c::gsl_test:: [15034]incorrect precision (%g obs vs %g expected)
-13, 1072693247, -350469331, 1058682594
FAIL: results.c::gsl_test:: [15036]incorrect precision (%g obs vs %g expected)
2142401953, 1061077358, 0, 0
FAIL: results.c::gsl_test:: [15038]incorrect precision (%g obs vs %g expected)
FAIL: results.c::gsl_test:: [15064]%s, %s (%g obs vs %g expected)
FAIL: results.c::gsl_test:: [15065]exceeded maximum number of iterations
FAIL: results.c::gsl_test:: [15066]incorrect precision (%g obs vs %g expected)
FAIL: results.c::gsl_test:: [15068]incorrect precision (%g obs vs %g expected)
FAIL: results.c::gsl_test:: [15070]incorrect precision (%g obs vs %g expected)
FAIL: results.c::gsl_test:: [15072]%s, %s
GSL roots test result is 1
 end of gsl roots test




GSL Examples
================
http://apwillis.staff.shef.ac.uk/aco/freesoftware.html  - passing function as argument
http://www.helsinki.fi/~fyl_tlpk/luento/ohj-13-GSL-e.html - passing function as argument
Fortran GSL
http://www.lrz.de/services/software/mathematik/gsl/fortran/index.html
http://uncwddas.googlecode.com/svn/trunk/other/gsl-1.8/gsl.vc/gsl.vc8.readme.txt


GSL Considerations
===================
32-64 bit
cblas replace with blas - can used ATLAS (C) at least


Build GSL with MinGW
========================
MSVC and MinGW Dlls - shows how to build Dll using MinGW that can be used by MSVC
http://www.mingw.org/wiki/MSVC_and_MinGW_DLLs

GSL can be build using MinGW with the following commands:
configure
OR  to just make a dynamic library without static lib;
./configure --enable-static=no
then
  make

In general, to specify on conditions and macro definitions:
env CPPFLAGS="-DGSL_DLL"  ./configure --enable-static=no
./configure --enable-static=no CPPFLAGS="-DGSL_DLL"
Don't need to add GSL_DLL definition - this is just an exmaple.


The result is are the following files:
in [gsl]\cblas\.libs:  libgslcblas.dll.a, libgslcblas.a, libgslcblas-0.dll, libgslcblas.la, libgslcblas.lai
in [gsl]\libs: libgsl.dll.a, libgsl.a, libgsl-0.dll, libgsl.la, libgsl.lai

If there are errors during MAKE, like these:
    infnan.c:98:3: error: #error "cannot define gsl_finite without HAVE_DECL_FINITE or HAVE_IEEE_COMPARISONS"
    infnan.c:115:3: error: #error "cannot define gsl_isnan without HAVE_DECL_ISNAN or HAVE_IEEE_COMPARISONS"
to solve this, go into [gsl]\config.h  and remove the UNDEF lines and add DEFINE lines like:
   //#undef HAVE_DECL_ISFINITE
   #define HAVE_DECL_ISFINITE 1

   //#undef HAVE_DECL_ISNAN
   #define HAVE_DECL_ISNAN 1


To make GSL compatible with MSVC(Microsoft Visual C), then edit gsl_types.h file, to ensure that the following line is defined:


 
Run-Time Check Failure #0 - The value of ESP was not properly saved
across a function call.  This is usually a result of calling a function
declared with one calling convention with a function pointer declared
with a different calling convention.

The message above about calling convention is misleading. The gsl dll may be in the right calling convention already. The caller to the functions in gsl need to specify explicitly
    __cdel(dllimport)
on the prototype for the function they need to call.


Running MSVC with minGW-built-GSL
===================================
First, built the GSL using MinGW as described in "Build GSL with MinGW" section.
Then,
To LINK with MSVC, use the *.dll.a, but rename them to *.lib
To RUN with MSVS, put the *-0.dll files into the runtime path.

Error Message when running
"Run-Time Check Failure #0 - The value of ESP was not properly saved across a function call.  This is usually a result of calling a function declared with one calling convention with a function pointer declared with a different calling convention."

..... need to check calling conventions.


Aside: To convert MSVC *.lib into *.a for GNU/GCC in MinGW, use the tool called "reimp".


GSL Makefile
=============
some key definitions of the Makefile is below.

$(SHELL) $(top_builddir)/libtool

libgsl.la: $(libgsl_la_OBJECTS) $(libgsl_la_DEPENDENCIES)
$(libgsl_la_LINK) -rpath $(libdir) $(libgsl_la_OBJECTS) $(libgsl_la_LIBADD) $(LIBS)


92:
libgsl_la_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
$(LIBTOOLFLAGS) --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \
$(libgsl_la_LDFLAGS) $(LDFLAGS) -o $@
294:
libdir = ${exec_prefix}/lib
343:
libgsl_la_LIBADD = $(SUBLIBS) $(am__append_1)
233:
LIBS = -lm
234:
LIBTOOL = $(SHELL) $(top_builddir)/libtool
112:
CCLD = $(CC)
184:
CFLAGS = -g -O2
344:
libgsl_la_LDFLAGS = -version-info $(GSL_LT_VERSION) $(am__append_2)


90:
am_libgsl_la_OBJECTS = version.lo
libgsl_la_OBJECTS = $(am_libgsl_la_OBJECTS)

88:
libgsl_la_DEPENDENCIES = $(SUBLIBS) $(am__append_1)
43:
am__append_1 = cblas/libgslcblas.la
44:
am__append_2 = -no-undefined
315:
SUBLIBS = block/libgslblock.la blas/libgslblas.la \
bspline/libgslbspline.la complex/libgslcomplex.la \
cheb/libgslcheb.la dht/libgsldht.la diff/libgsldiff.la \
deriv/libgslderiv.la eigen/libgsleigen.la err/libgslerr.la \
fft/libgslfft.la fit/libgslfit.la histogram/libgslhistogram.la \
ieee-utils/libgslieeeutils.la integration/libgslintegration.la \
interpolation/libgslinterpolation.la linalg/libgsllinalg.la \
matrix/libgslmatrix.la min/libgslmin.la monte/libgslmonte.la \
multifit/libgslmultifit.la multimin/libgslmultimin.la \
multiroots/libgslmultiroots.la ntuple/libgslntuple.la \
ode-initval/libgslodeiv.la ode-initval2/libgslodeiv2.la \
permutation/libgslpermutation.la \
combination/libgslcombination.la multiset/libgslmultiset.la \
poly/libgslpoly.la qrng/libgslqrng.la randist/libgslrandist.la \
rng/libgslrng.la roots/libgslroots.la siman/libgslsiman.la \
sort/libgslsort.la specfunc/libgslspecfunc.la \
statistics/libgslstatistics.la sum/libgslsum.la \
sys/libgslsys.la test/libgsltest.la utils/libutils.la \
vector/libgslvector.la cdf/libgslcdf.la \
wavelet/libgslwavelet.la


Sample of actual output in Make:

libtool: link: ln .libs/libgsl.lax/libgslcdf.a/fdist.o .libs/libgsl.lax/lt72-fdist.o || cp .libs/libgsl.lax/libgslcdf.a/fdist.o .libs/libgsl.lax/lt72-fdist.o
libtool: link: ln .libs/libgsl.lax/libgslcdf.a/flat.o .libs/libgsl.lax/lt73-flat.o || cp .libs/libgsl.lax/libgslcdf.a/flat.o .libs/libgsl.lax/lt73-flat.o
libtool: link: ln .libs/libgsl.lax/libgslcdf.a/gamma.o .libs/libgsl.lax/lt74-gamma.o || cp .libs/libgsl.lax/libgslcdf.a/gamma.o .libs/libgsl.lax/lt74-gamma.o
libtool: link: ln .libs/libgsl.lax/libgslcdf.a/gauss.o .libs/libgsl.lax/lt75-gauss.o || cp .libs/libgsl.lax/libgslcdf.a/gauss.o .libs/libgsl.lax/lt75-gauss.o
libtool: link: ln .libs/libgsl.lax/libgslcdf.a/geometric.o .libs/libgsl.lax/lt76-geometric.o || cp .libs/libgsl.lax/libgslcdf.a/geometric.o .libs/libgsl.lax/lt76-geometric.o
libtool: link: ln .libs/libgsl.lax/libgslcdf.a/laplace.o .libs/libgsl.lax/lt77-laplace.o || cp .libs/libgsl.lax/libgslcdf.a/laplace.o .libs/libgsl.lax/lt77-laplace.o
libtool: link: ln .libs/libgsl.lax/libgslcdf.a/logistic.o .libs/libgsl.lax/lt78-logistic.o || cp .libs/libgsl.lax/libgslcdf.a/logistic.o .libs/libgsl.lax/lt78-logistic.o



General C define macros and typedef
=====================================
#define [dummy] [real text]
eg.
# define __BEGIN_DECLS extern "C" {

this means whenever the "__BEGIN_DECLS" is used in the code, the pre-processor will replace it with:
extern "C" {


typedef [type] [alias]
eg.

/* type declaration */
typedef struct {
int number;
char *text;
} LINE;


/* Variable declaration */
LINE buffer[MAXLINES];


Note the actual type, like the struct above, is declared first, followed by the last word which is the alias LINE.
BUT, this is reverse for #define where the actual alias or dummy is the first word, then the rest are the real text.

An alternative way to declare a struct using "tag" is:

/* type declaration */
struct templ {
int number;
char *text;
} ;
/* Variable declaration */
struct templ buffer[MAXLINES];


Note that using the tag method for struct, the keyword struct still need to be included in the declaration.



GSL-C-Fortran Framework
=======================
1. Compile the GSL into a library with vcpp_gsl.lib|.bin, see "Building using Visual Studio C/C++"
2. Create a wrapper C++ project, say gslCwrap.
3. Create a header file gslCwrap.h and code file gslCwrap.c
4. In the header file gslCwrap.h
a) prepare to export the wrappers.

extern "C"{
__declspec(dllexport) void cdfGaussP (double x, double *p);
__declspec(dllexport) void fnGamma (double x, double *gamma);
}

b) include reference to GSL header like and call gsl function:

#include "gsl/gsl_cdf.h"
#include "gsl/gsl_sf_gamma.h"
void cdfGaussP (double x, double *p) {      *p =  gsl_cdf_ugaussian_P(x); }// 2 * x;  }
void fnGamma (double x, double *gamma) {    *gamma = gsl_sf_gamma(x); }

the actual gsl function names used here are: gsl_cdf_ugaussian_P, gsl_sf_gamma


5. To call from Fortran, prepare a Fortran interface:

    interface
        SUBROUTINE cdfGauss_P(x, P)  BIND(C, NAME='cdfGaussP')
            USE, INTRINSIC :: ISO_C_BINDING     ! New Fortran 2003 standard
            IMPLICIT NONE
            REAL(C_DOUBLE), VALUE :: x
            REAL(C_DOUBLE) :: P
        END SUBROUTINE cdfGauss_P
        SUBROUTINE fnGamma(x, dGamma)  BIND(C, NAME='fnGamma')
            USE, INTRINSIC :: ISO_C_BINDING     ! New Fortran 2003 standard
            IMPLICIT NONE
            REAL(C_DOUBLE), VALUE :: x
            REAL(C_DOUBLE) :: dGamma
        END SUBROUTINE fnGamma
end interface

6. In the Fortran code, just call the function as specified by the interface.
Eg.
    call fnGamma(vecX(ii), answ)      
call cdfGauss_P(vecQ(ii), answ)      


There are at least 3 projects involved here.
vcpp_gsl - this keeps the gsl code in its own project, so that it is independent of any change in interface or usage.
gslCwrap - this consist of the C wrappers so that it can be used by C code or Fortran code. It also transforms all gsl functions into VOID C functions so that it can be called as Fortran subroutines. Fortran Interface - this is freely linked to the gslCwrap wrapper code. There is no direct reference to the GSL code.  

Tuesday, January 11, 2011

Callback to C# from Unmanaged Fortran - PASSING ARRAYS

Thanks to reader kensun87 for requesting this feature.

This article builds on another article which should be read first:
Callback to C# from Unmanaged Fortran
http://xtechnotes.blogspot.com/2008/07/callback-to-c-from-unmanaged-fortran.html


The original article on callback referenced above passes a single value (scalar) between callbacks. This example passes an ARRAY of integers. This is much trickier because arrays need to be handled using IntPtr when passing between callbacks.

The Fortran code and C# code is listed below first, then some explanations. Please note that the fundamental explanations on callback will not be here, instead please see the previous article referenced above.

-------  Fortran code  ------

module f90Callback
    contains


    ! used by CScallbackDriver::Program.cs    
    subroutine func1(iArr, progressCllBak)
    !DEC$ ATTRIBUTES DLLEXPORT ::func1
    !DEC$ ATTRIBUTES REFERENCE :: iArr, progressCllBak
        implicit none
        external progressCllBak
        integer :: iCB
        integer, INTENT(OUT) :: iArr(2)
        
     ! Body of f90Callback
        print *, "Hello from Fortran func1(iArr, progressCllBak)"  
        iCB = 3
        iArr(1) = 5
        iArr(2) = 7
        call progressCllBak(iArr, 2)
        print *, "Final Fortran arrays are ", iArr
         
        return
    end subroutine func1
    
end module f90Callback


-------  C# code  ------

using System;
using System.Collections.Generic;
using System.Text;
using System.Runtime.InteropServices;


namespace CScallbackDriver
{
    unsafe class Program
    {
        // 0. Define a counter for the Progress callback to update
        public int localCounter; 


        // 1. Define delegate type
        [UnmanagedFunctionPointer(CallingConvention.Cdecl)]
        public delegate void dgateIntPtr( IntPtr numYears, ref int iSize);


        // 2. Create a delegate variable
        public dgateIntPtr dg_progCBPtr;


        public Program() {
            // 3. Instantiate delegate, typically in a Constructor of the class
            dg_progCBPtr = new dgateIntPtr(onUpdateProgressPtr);
        }


        // 4. Define the c# callback function
        ///

        /// Callback function where an array is passed from Fortran
        ///

        /// pointer to the unmanaged array
        /// size of the array represented by progCount
        public void onUpdateProgressPtr( IntPtr progCount, ref int iSize)
        {
            //Unsafe code
            unsafe
            {
                Console.WriteLine("In C# callback function");


// reading in array from unmanaged code
                int[] managedArr2 = new int[2];
                Marshal.Copy(progCount, managedArr2, 0, iSize);
                
// writing out array to unmanaged code
                managedArr2[1] = 99;
                Marshal.Copy(managedArr2, 0, progCount, iSize);
                Console.WriteLine("Going out of C# callback function");
            }


        }
         
        static void Main(string[] args)
        {
            Program myProg = new Program();   
            myProg.localCounter = 0;
            int[] iArrB = new int[2];
            
            //6. Call normal Fortran function from DLL, and passing the callback delegate
            func1Ptr(ref iArrB[0], myProg.dg_progCBPtr);


            Console.ReadKey();   
        }


        // 5. Define the dll interface 
        // Pointer
        [DllImport("f90Callback", EntryPoint = "F90CALLBACK_mp_FUNC1", CharSet = CharSet.Auto, CallingConvention = CallingConvention.Cdecl)]
        public static extern void func1Ptr([In, Out] ref int iArr, [MarshalAs(UnmanagedType.FunctionPtr)]  dgateIntPtr blah);
    }
}






A few things to note:
1. in the Main method, you can ignore iArrB as it does not have implications in the callback.
2. In the Fortran code, when the call to the C# callback is made, note it is important to pass the correct size of the array, in this example 2:
              call progressCllBak(iArr, 2)
3. In the Fortran code, the iArr is initialised with values {5,7}
4. In C# code, in step 1, see exactly how the arguments of the delegate dgateIntPtr are defined.
5. In C# code, in step 4, this is where the actual callback function is defined.
- the name of the callback function is onUpdateProgressPtr with arguments: ( IntPtr progCount, ref int iSize)
- the array called managedArr2 is defined with size iSize.
- Marshal.Copy is used to copy the unmanaged array progCount, into C# array managedArr2.
- the second element of the array is changed to 99.
- Marshal.Copy is used again but to copy the opposite way from managed array manageArr2 to unmanaged array progCount.
6. The call back returns to the Fortran function which then prints out the array where second element has new value of 99.
7. When the Fortran function finishes, it returns control back to C# main method.
8. All other steps in this example are similar to the previous callback example where only a scalar integer is passed.

Thursday, November 18, 2010

Matlab Function in Fortran - perms - permutation

Below is the Fortran implementation of the perms algorithm found in Matlab. It creates a matrix whose rows consist of all possible permutations of the n elements of vector vecIn.

There are many useful utility functions in Matlab and occassionally an important but rare (difficult to find source in the internet) algorithm such as this permutation algorithm. Hopefully there will be more Fortran version of Matlab utility functions appearing here in the future.

The one listed below is unpolished and unoptimized and not suitable for large inputs. It is a very basic implementation taken directly from the mathematical definition. It is also tested only for a small number of cases. Corrections, suggestions and comments are welcomed.




!/*
! * perms - All Possible permutations
! * 
! *     (C)  Copyright 2010 xTechNotes.blogspot.com
! * 
! *   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 3 of the License, or
! *   (at your option) any later version.
! *
! *    This program is distributed in the hope that it will be useful,
! *    but WITHOUT ANY WARRANTY; without even the implied warranty of
! *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! *    GNU General Public License for more details.
! *
! *    You should have received a copy of the GNU General Public License
! *    along with this program.  If not, see .
! */ 


    RECURSIVE subroutine perms(n, nfact, vecIn, vecOut)
        IMPLICIT NONE 
        INTEGER, INTENT(IN) :: n                ! number of elements of original vector
        INTEGER, INTENT(IN) :: nfact            ! Factorial(n)
        REAL, INTENT(IN) :: vecIn(n)            ! original vector to be permuted
        REAL, INTENT(OUT):: vecOut(nfact, n)    ! Permutations of original vector


        ! local variables         
        real :: vecInTmp(n)
        integer :: ii , mfact
        real :: dtmp
        
        vecOut = 0.0
        if (n .le. 0) then 
            return
        elseif ( n .eq. 1 ) then            
            vecOut(1,1) = vecIn(1)
        elseif ( n .eq. 2 ) then
            vecOut(1, :) = vecIn(:)
            vecOut(2, :) = (/vecIn(2), vecIn(1)/)
        else        
            ! ii = 1
            call factorial(n-1, mfact)
            vecOut(1:mfact, 1) = vecIn(1)
            call perms(n-1, mfact,  vecIn(2:n), vecOut(1:mfact, 2:n))
            
            do ii = 2, n
                vecInTmp = VecIn
                vecInTmp(1) = VecIn(ii)
                vecInTmp(ii) = VecIn(1)
                vecOut((ii-1)*mfact+1 : ii*mfact, 1) = vecInTmp(1)                
                call perms(n-1, mfact, vecInTmp(2:n), vecOut( (ii-1)*mfact+1 : ii*mfact, 2:n))
            enddo 
        
        endif    
        
    end subroutine perms

Friday, November 12, 2010

Matlab Function in Fortran - conv2 - Convolution in 2D

Below is the Fortran implementation of the Convolution 2D algorithm, or conv2 as found in Matlab. There are many useful utility functions in Matlab and occassionally an important but rare (difficult to find source in the internet) algorithm such as this convolution algorithm. Hopefully there will be more Fortran version of Matlab utility functions appearing here in the future.

The one listed below is unpolished and unoptimized and not suitable for large inputs. It is a very basic implementation taken directly from the mathematical definition of discrete convolution in 2D. It is also tested only for a small number of cases. Corrections, suggestions and comments are welcomed.




!/*
! * CONV2 - Convolution in 2D
! *
! *     (C)  Copyright 2010 xTechNotes.blogspot.com
! *
! *   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 3 of the License, or
! *   (at your option) any later version.
! *
! *    This program is distributed in the hope that it will be useful,
! *    but WITHOUT ANY WARRANTY; without even the implied warranty of
! *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! *    GNU General Public License for more details.
! *
! *    You should have received a copy of the GNU General Public License
! *    along with this program.  If not, see .
! */  
  
    subroutine conv2(nXi, nXj, nHi, nHj, nYi, nYj, X, H, Y)
        INTEGER, INTENT(IN) :: nXi, nXj, nHi, nHj           ! Dimensions of Input, Kernel
        INTEGER, INTENT(IN) :: nYi, nYj                     ! Dimensions of Output - MUST BE nXi+nHi-1, nXj+nHj-1
        REAL, INTENT(IN) :: X(0:nXi-1, 0:nXj-1)   ! Input
        REAL, INTENT(IN) :: H(0:nHi-1, 0:nHj-1)   ! Kernel
        REAL, INTENT(OUT) :: Y(0:nYi-1, 0:nYj-1)  ! Output
        integer :: im, in, ii, ij, ip, iq
      
        Y = 0.0d0
        if( nYi .ne. nXi+nHi-1  .or. nYj .ne. nXj+nHj-1 ) RETURN
      
        do in = 0, nYj-1
            do im = 0, nYi-1
                do ij = 0, nXj-1
                    do ii = 0, nXi-1
                    ip = im - ii
                    iq = in - ij
                    if (ip .ge. 0 .and. ip .le. nHi-1 .and. iq .ge. 0 .and. iq .le. nHj-1) then
                        Y(im, in) = Y(im, in) + X(ii, ij) * H(ip, iq)
                    endif                      
                    enddo
                enddo
            enddo
        enddo
        return
    end subroutine conv2

Thursday, June 24, 2010

Fortran Debugging, Threading, Optimising articles

This post is a quick summary to some online articles which I find useful. Only the extract are given below. The full articles can be found in their original sources via their URL.


Threading Fortran applications for parallel performance on multi-core systems

Most processors now come with multiple cores, and future increases in performance are expected to come mostly from increases in core count. Performance sensitive applications that neglect the opportunities presented by additional cores will soon be left behind. This article discusses ways for an existing, serial Fortran application to take advantage of these opportunities on a single, shared memory system with multiple cores. Issues addressed include data layout, thread safety, performance and debugging. Intel provides software tools that can help in the development of robust, parallel applications that scale well.

Levels of Parallelism

1 SIMD instructions
2 Instruction level
3 Threading (usually shared memory)
4 Distributed memory clusters
5 “embarassingly parallel” multiprocessing

Ways to Introduce Threading

1 Threaded libraries, e.g. Intel® MKL
2 Auto-parallelization by the compiler
3 Asynchronous I/O (very specialized; see compiler documentation)
4 Native threads
5 OpenMP

Intel® Math Kernel Library
1 Many components of MKL have threaded versions
2 Link threaded or non-threaded interface
3 Set the number of threads

Example: PARDISO (Parallel Direct Sparse Solver)
1 Solver for large, sparse symmetric and antisymmetric systems of linear equations on shared memory systems
2 For algorithms, see http://www.pardiso-project.org

Auto-parallelization

1 The compiler can thread simple loops automatically
2 Based on the same RTL threading calls as OpenMP:

Conditions for Auto-parallelization

1 Loop count known at entry (no DO WHILE)
2 Loop iterations are independent
3 Enough work to amortize parallel overhead
4 Conditions for OpenMP loops are similar
5 Directives may be used to guide the compiler:

Example: matrix multiply


OPENMP - advantages

1 Standardized API based on compiler directives

OpenMP Programming Model

Fork-Join Parallelism:

1 Master thread spawns a team of threads as needed

Note that Intel’s implementation of OpenMP creates a separate monitor thread in addition to any user threads.


OPENMP – where to thread

1 Start by mapping out high level structure
2 Where does your program spend the most time?
3 Prefer data parallelism
4 Favor coarse grain (high level) parallelism

Example: Square_Charge

1 calculates the electrostatic potential at a series of points in a plane
   due to a uniform square distribution of charge


Openmp: How do threads interact?

1 OpenMP is a shared memory model
2 Unintended sharing of data causes race conditions:
3 To control race conditions…
4 Synchronization is expensive so…

OPENMP – data

1 Identify which data are shared between threads, which need a separate copy for each thread

2 It’s helpful (but not required) to make shared data explicitly global in Modules or common blocks,
   thread private data as local and automatic.

3 Dynamic allocation is OK (malloc, ALLOCATE)

4 Each thread gets its own private stack, but the heap is shared by all threads

OPENMP – data scoping

1 Distinguish lexically explicit parallel regions from the “dynamic extent” (Functions or subroutines called from within an explicit parallel region. These might contain no OpenMP directives or only “orphaned” OpenMP directives)

2 Lexically explicit: !$OMP PARALLEL to !$OMP END PARALLEL


Thread Safety

1 A threadsafe function can be called simultaneously from multiple threads, and still give correct results
2 ifort serial defaults:
3 When compiling with –openmp, default changes

Making a function thread safe

1 With the compiler
2 In source code
3 In either case:
4 OpenMP has various synchronization constructs to protect operations that are potentially unsafe

Thread Safe Libraries

1 The Intel® Math Kernel library is threadsafe
2 The Intel Fortran runtime library has two versions

Performance considerations

1 Start with optimized serial code, vectorized inner loops, etc. (-O3 –xsse4.2 –ipo …)
2 Ensure sufficient parallel work
3 Minimize data sharing between threads
4 Avoid false sharing of cache lines
5 Scheduling options

Timers for threaded apps

1 The Fortran standard timer CPU_TIME returns “processor time”
2 The Fortran intrinsic subroutine SYSTEM_CLOCK returns data from the real time clock
3 dclock (Intel-specific function) can also be used

Thread Affinity Interface

1 Allows OpenMP threads to be bound to physical or logical cores

NUMA considerations

1 Want memory allocated “close” to where it will be used
2 Remember to set KMP_AFFINITY


Common problems

1 Insufficient stack size
2 For whole program (shared + local data):
3 For individual thread (thread local data only)

Tips for Debugging OpenMP apps

1 Run with OMP_NUM_THREADS=1
2 Build with -openmp-stubs -auto
3 If works without –auto, implicates changed memory model
4 If debugging with PRINT statements
5 Debug with –O0 –openmp

Floating-Point Reproducibility

1 Runs of the same executable with different numbers of threads may give slightly different answers
2 Floating point reductions are still not strictly reproducible in OpenMP, even for same number of threads

Intel-specific Environment Variables

1 KMP_SETTINGS = 0 | 1
2 KMP_VERSION = off | on
3 KMP_LIBRARY = turnaround | throughput | serial
4 KMP_BLOCKTIME
5 KMP_AFFINITY (See main documentation for full API)
6 KMP_MONITOR_STACKSIZE
7 KMP_CPUINFO_FILE

Tools for Debugging OpenMP apps

1 The compiler source checker (‘parallel lint’)
2 Updated Intel Parallel Debugger, idb (Linux) and Intel Parallel Debugger Extension (on Windows)

Intel® Thread Checker
1 Unified set of tools that pinpoint hard-to-find errors in multi-threaded applications
2 Display data at the Linux command line or via a Windows GUI

Intel® Thread Profiler
1 Features & Benefits

Summary
Intel software tools provide extensive support for threading applications to take advantage of multi-core architectures.
Advice and background information are provided for a variety of issues that may arise when threading a Fortran application.





Tips for Debugging Run-time Failures in Applications Built with the Intel(R) Fortran Compiler

Your app builds successfully, but crashes at runtime. What Next? Try some useful Intel compiler diagnostic options before launching into lengthy debugger sessions.

1) Build with /traceback (Windows*) or –traceback (Linux* or Mac OS* X).

2) Build with /gen-interfaces /warn:interfaces (Windows) or –gen-interfaces –warn interfaces (Linux or Mac OS X).

3) Try building and running with /check (Windows) or –check (Linux and Mac OS X).

4) Build your program, including the main routine, with /fpe:0 (Windows) or –fpe0 (Linux or Mac OS X).

5) If your application fails early on with a segmentation fault, you might be exceeding the default maximum stack size. On Linux or Mac OS X, try setting
ulimit –s unlimited (bash) or limit stacksize unlimited (C shell)

6) Use the compiler provided interfaces. If you call run-time library functions, build with
      USE IFLPORT
If you call OpenMP run-time library functions, compile with
      USE OMP_LIB
If you call functions from MKL or IMSL*, USE the corresponding module(s).

7) Look carefully for any error messages in your output log file.

8) If you are building an application using OpenMP*, check out the advice under “Tips for Debugging OpenMP Apps” at http://software.intel.com/en-us/articles/threading-fortran-applications-for-parallel-performance-on-multi-core-systems/


9) For Windows, see the section Building Applications / Debugging in the main compiler documentation. For Linux or Mac OS X, see the documentation for the Intel(R) Debugger (idb).

Saturday, November 14, 2009

Notes 64bit

Notes64bit

Contents
============
References
Definition
Article - The 64-Bit Advantage
Article - x86: registered offender
Itanium2
Feature Comparison
Registers
AMD K8 vs Conroe FPU
Porting to a 64-bit Intel® architecture
How to check if code is 32bit or 64bit
Limitations
Large Arrays



References
=============
http://www.intel.com/cd/ids/developer/asmo-na/eng/197664.htm?page=5


Definition
===========
Ref: http://en.wikipedia.org/wiki/64-bit
"64-bit" computer architecture generally has integer registers that are 64 bits wide, which allows it to support (both internally and externally) 64-bit "chunks" of integer data.



Size to consider are: registers, address buses, or data buses.


Most modern CPUs such as the Pentium and PowerPC have 128-bit vector registers used to store several smaller numbers, such as 4 32-bit floating-point numbers. A single instruction can operate on all these values in parallel (SIMD). They are 128-bit processors in the sense that they have 128-bit registers and in some cases a 128-bit ALU, but they do not operate on individual numbers that are 128 binary digits in length.


Article - The 64-Bit Advantage
===============================
Ref: http://www.pcmag.com/print_article/0,3048,a=116259,00.asp

The 32-bit Pentium-class chips that dominate today's desktops fetch and execute instructions from system memory in 32-bit chunks; 64-bit chips handle 64-bit instructions. And that's just what the workstation-class Intel Itanium 2 and HP Alpha chips do inside the TeraGrid's clusters.

New desktop-class 64-bit chips, such as the AMD Athlon64 and the Apple/ IBM PowerPC G5, can handle 64-bit instructions as well, but most PC apps—even the few that optimize some operations to exploit 64-bit processing—still rely on 32-bit instructions. A new generation of games and apps will no doubt take fuller advantage of 64-bit chips. But their ability to harness the new architecture fully may be hampered by the need to interact with Windows, since none of the desktop versions of the OS is yet slated for 64-bit optimization.

A major advantage to 64-bit processors over their 32-bit cousins is support for greater amounts of memory. In theory, a 64-bit processor can address exabytes (billions of billions of bytes) of RAM; 32-bit chips can use a maximum of 8GB of RAM. This breakthrough is used to good advantage at the National Center for Supercomputing Applications' (NCSA) TeraGrid, which allocates 12GB of system memory each to half of its 256 Itanium 2 processor nodes. It will be a while before anyone knows how fast Quake would run with that much memory, since PC motherboards don't exceed 8GB of RAM.

Future 64-bit apps will be able to chew on a class of computations known as floating-point operations far faster than 32-bit apps can. Necessary for 3-D rendering and animation of everything from molecular models to Halo aliens, floating-point calculations are so essential to complex scientific analysis that FLOPS (floating-point operations per second) are used as the unit of supercomputing performance. The ability of 64-bit chips to process floating-point operations faster and far more precisely than their 32-bit counterparts make them powerhouses for simulations and visualization.


Article - x86: registered offender
===================================
Ref: http://techreport.com/reviews/2005q1/64-bits/index.x?pg=2

"Another problem with the x86 ISA is the number of general-purpose registers (GPRs) available. Registers are fast, local slots inside a processor where programs can store values. Data stored in registers is quickly accessible for reuse, and registers are even faster than on-chip cache. The x86 ISA only provides eight general-purpose registers, and thus is generally considered register-poor. Most reasonably contemporary ISAs offer more. The PowerPC 604 RISC architecture, to give one example, has 32 general-purpose registers. Without a sufficient number of registers for the task at hand, x86 compilers must sometimes direct programs to spend time shuffling data around in order to make the right data available for an operation. This creates overhead that slows down computation.

To help alleviate this bottleneck, the x86-64 ISA brings more and better registers to the table. x86-64 packs 8 more general-purpose registers, for a total of 16, and they are no longer limited to 32-bit values—all 16 can store 64-bit datatypes. In addition to the new GPRs, x86-64 also includes 8 new 128-bit SSE/SSE2 registers, for a total of 16 of those. These additional registers bring x86 processors up to snuff with the competition, and they will quite likely bring the largest performance gains of any aspect of the move to the x86-64 ISA.

What is the magnitude of those performance gains? Well, it depends. Some tasks aren't constrained by the number of registers available now, while others will benefit greatly when recompiled for x86-64 because the compiler will have more slots for local data storage. The amount of "register pressure" presented by a program depends on its nature, as this paper on 64-bit technical computing with Fortran explains:

The performance gains from having 16 GPRs available will vary depending on the complexity of your code. Compute-intensive applications with deeply nested loops, as in most Fortran codes, will experience higher levels of register pressure than simpler algorithms that follow a mostly linear execution path. "

Summary -
x86 - 8x 32-bit General Purpose Registers
x86-64 - 16x 64-bit General Purpose Registers
Fortran - more do loops need bigger and more GPRs


Itanium2
=========
Ref: http://www.itmanagersjournal.com/feature/8611

Intel's Itanium 2, or IA64, is unlike any of the other 64-bit processors in production. It uses a Very Long Instruction Word (VLIW) design that depends on the software's compiler for performance. When the compiler creates program binaries for the Itanium 2, it predicts the most efficient method of execution, so the processor does less work when the program is running -- the software schedules its own resources beforehand, rather than forcing the hardware to do it on the fly. IA64 is used in the same kinds of workstations that UltraSPARC processors are used in, and can also scale up to 128 processors in high-powered servers. Silicon Graphics and Hewlett-Packard both sell computers based on the Itanium 2. GNU/Linux is generally the operating system of choice for IA64-based systems, but HP-UX and Windows 2003 Server will work on HP Itanium 2 servers.



Feature Comparison
==================

Size of Fetch and Execute Instructions - 64bit vs 32bit chunks
Number of General Purpose Registers (Fetch Registers?)
Memory Access - 18.4x10^9 GB vs 4GB
Floating Point Operations - faster with 64bit than 32bit

Vector Registers - eg Pentium has 128bit data registers which store up to 4 32bit data register.
ALU
FPU


Itanium2
- good for FP processing
- 2FPU (=1 FMAC or 2multiplication and 1 add), plus additional 2 FMACs for 3D processing.
- 64 bit address space
- a derivative of VLIW, dubbed Explicitly Parallel Instruction Computing (EPIC). It is theoretically capable of performing roughly 8 times more work per clock cycle than a non-superscalar CISC or RISC architecture due to its Parallel Computing Microarchitecture.
- support 128 integer, 128 floating point, 8 branch and 64 predicate registers (for comparison, IA-32 processors support 8 registers and other RISC processors support 32 registers


UltraSparc T1
- 8x Integer Cores share 1 FPU
- good for integer processing compared to Itanium


Throughout its history, Itanium has had the best floating point performance relative to fixed-point performance of any general-purpose microprocessor. This capability is not needed for most enterprise server workloads. Sun's latest server-class microprocessor, the UltraSPARC T1 acknowledges this explicitly, with performance dramatically skewed toward the improvement of integer processing at the expense of floating point performance (eight integer cores share a single FPU). Thus Itanium and Sun appear to be addressing separate subsets of the market. By contrast, IBM's cell microprocessor, with a single general-purpose POWER core controlling eight simpler cores optimized for floating point, may eventually compete against Itanium for floating-point workloads.


Registers
==========
Integer - can be used to store pointers
Floating Point - most CPUs also have FPUs
Other

examples:
x86 - has x87 FPU with 8 x 80bit registers
x86 with SSE - 8x 128bit FP registers
x86-64 - has SSE with 16x 128bit FP registers
Alpha - has 32x 64bit FP registers and 32x 64bit integer registers.
Itanium2 - 128x 64bit GPRs, 128x 82bit FPregisters, 64x 1bit predicates, 8x 64bit branch registers

AMD K8 vs Conroe FPU
======================
Possibly better floating point performance of K8 processors

http://www.xbitlabs.com/articles/cpu/display/amd-k8l_5.html


Porting to a 64-bit Intel® architecture
============================================
(ref: http://www.developers.net/intelisnshowcase/view/358)

Porting Application Source Code
The most significant issues that software developers should face in porting source code to the 64-bit world concern the changes in pointer size and fundamental integer types. As such, these differences should appear most prominently in C and C++ programs. Code written in Fortran, COBOL, Visual Basic, and most other languages (except assembly language, which must be completely rewritten), will need no modification. A simple recompilation is often all that is needed. Java code should not even need recompilation; Java classes should execute the same on a 64-bit JVM as on any 32-bit virtual machine.

C (from here on, C++ is included in all discussions of C) code, however, by allowing casting across types and direct access to machine-specific integral types will need some attention.

The first aspect is the size of pointers; 64-bit operating systems use 64-bit pointers. This means that the following will equal eight (and no longer four):

sizeof (ptrdiff_t)

As a result, structures that contain pointers will have different sizes as well. As such, if data laid out for these structures is stored on disk, reading it in or writing it out will cause errors. Likewise, unions with pointer fields will have different sizes and can cause unpredictable results.

The greatest effect, though, is felt wherever pointers are cast to integral types. This practice, which has been condemned for years as inimical to portability, will come back to haunt programmers who did not abandon it. The problems caused by it are traceable to the different widths used by pointers, integers and longs on the various platforms. Let's examine these.


How to check if code is 32bit or 64bit
========================================
use the dumpbin utility and look for the output under FILE HEADER VALUES.
eg.
dumpbin /headers hello.exe

Results if 64 bit:
FILE HEADER VALUES
8664 machine (x64)

Results if 32 bit:
FILE HEADER VALUES
14C machine (x86)

Limitations
=============
Virtual Address Limit - theoretical 16EB
Virtual Address Limit - practical
i) Windows use 44bits -> 16TB, apparently allow only 8TB to be used.


Large Arrays
=============
http://episteme.arstechnica.com/eve/forums/a/tpc/f/6330927813/m/420003239831/r/308002539831

BigArray, getting around the 2GB array size limit
http://blogs.msdn.com/joshwil/archive/2005/08/10/450202.aspx

Monday, January 05, 2009

How to use Brook+ for GPU computing

AMD Stream Computing allows developers to use the GPU to perform parallel computations for HPC applications. This guide is meant to complement the AMD Stream Computing User Guide. It is essential to read the official User Guide to gain a brief understanding before following the notes below.

Ref: http://ati.amd.com/technology/streamcomputing/Stream_Computing_User_Guide.pdf

System - The following notes are compiled based on the following system.
Intel CPU
ATI Radeon (Check cards for GPU computing capability)
Microsoft Visual Studio .Net with C/C++ compilers - for C/C++ code
Intel Visual Fortran Compilers - for Fortran code
Brook+ SDK by AMD - to compile Brook code


br source file
================
The Brook+ source file contain code that follow C/C++ syntax and is compiled/pre-processed by the Brook+ compiler into C/C++ file. Both Brook+ functions and C/C++ functions can exist within the same *.br file. The Brook+ functions are the functions that utilises the GPU hardware.

Example of a Brook+ function is given below:
kernel void sumaa(float a<>, float b<>, out float c<>){
c = a + b;
}

1. Special Brook+ keywords (ie. Not C/C++ words): kernel, out
2. Note the template like structures "float a<>" which are recognized by the Brook+ compiler. They indicate stream / GPU data type and are not the same as C++ templates.
3. Multiple functions like the above can exist in the same *.br file. Other normal C/C++ functions can also exist inside the *.br file.


Compiling Brook+ Code (*.br)
=============================
0. Open up a Command Console and go to the directory where the *.br file is located.

1. To compile code called sum.br:
\sdk\bin\brcc_d -k sum.br
where is the installation directory of the Brook SDK from AMD.

2. This the brook+ compiler / preprocessor creates the following in the same directory.
sum.cpp
sum.h
sum_gpu.h

3. A few notes to consider
i) There are two compilers: brcc and brcc_d. They correspond to brook.lib/dll and brook_d.lib/dll respectively.
Using the wrong combination may crash the program during execution.
ii) The -k option generates intermediate code that may be useful for use with the AMD's Stream Kernel Analyzer.
iii) The C/C++ code that are generated need to be compiled using standard C/C++ compilers and link to the proper libraries and dlls, hence the next section.
iv) Before v1.3, C/C++ wrapper functions, also known as host side code, exist within the *.br source file. As of v1.3, the host side code can be written in C++ and exist in a separate normal C++ file, provided it is configured with the proper include and lib directory information.



Compiling the C/C++ code
==========================
This step produces a win32 DLL from the C/C++ code that are generated by Brook+. The resultant DLL should be
able to be used by other win32 applications (eg C++ or Fortran).

1. From Visual Studio .Net, Open a new solution / project by:
Add Project -> Visual C++ -> Win32 -> Win32 project.
In the Application Settings dialog, select DLL, Export Symbols

2. Add the *.br and the files generated by the Brook+ compiler into the current project by using
"Add existing file".

3. Under the Project Property configuration pages, add the following settings:
C++ -> Additional Include Directories: \sdk\include
C++ -> Code Generation -> Runtime Library: Multi-threaded Debug DLL (/MDd)
C++ -> Advanced -> Calling Convention: __cdecl (/Gd)
Linker -> Additional Library Directories: \sdk\lib
Linker -> Input -> Additional Dependencies: \sdk\lib\brook_d.lib

4. When the *.br is modified, compile the *.br files in Command Console, then compile the generated c/c++ code from within the VisualStudio.Net environment.

Some Notes:
i) One can configure VisualStudio.Net to accept *.br files and compile using the Brook+ compiler. However, I find
that it still requires the user to manually initiate compilation for Brook files and then for C/C++ files. Hence,
I don't find it to be any efficient than compiling by command line.
ii) The *.br source files can be added to the project and can be edited using the VS.Net environment.


The C/C++ driver or library wrapper
====================================
The Brook+ functions need to be wrapped or called directly from C/C++ functions. For the purpose of creating DLL functions, we will put C/C++ wrappers over the Brook+ functions.

The usage of the Brook+ functions involve 3 steps. Each of these step are described with examples here:

Declaring and sizing variables - the meaning and reason for the declarations will become clear in the following sections.
// Normal C/C++ variables
float input_a[10][10];
float input_b[10][10];
float input_c[10][10];
float input_a1[10];
float input_b1[10];
float input_c1[10];

// For dimensioning Brook+ variables
unsigned int ileng = 10;
unsigned int dims[2] = {10,10};
unsigned int dim1[1] = {10};

// Equivalent Brook+ variables
brook::Stream a(2, dims);
brook::Stream b(2, dims);
brook::Stream a1(1, dim1);
brook::Stream b1(1, dim1);
brook::Stream *c = new brook::Stream(1, &ileng);
brook::Stream *d = new brook::Stream(2, dims);
brook::Stream d1(1, dim1);


// Assign values to normal C/C++ vectors and matrices for:
// input_a1, input_b1, input_a, input_b
..................

1. Reading normal C/C++ variables into Brook+ variables
a.read(input_a);
b.read(input_b);
a1.read(input_a1);
b1.read(input_b1);
This step transforms a normal C/C++ variable into a Brook+ variable which the GPU can understand. No other manipulation need to be done to the Brook+ variable.

2. Performing the computation by calling the Brook+ function
sumaa(a,b,*d); // operating on a matrix
sumaa(a1,b1,d1); // operating on a vector

3. Writing the output from Brook+ into normal C/C++ variables
// old method
streamWrite(*d, input_c);
streamWrite(d1, input_c1);
// new method
d->write(input_c);
d1.write(input_c1);
Once the Brook+ variable has been copied back to a normal C/C++ variable, one can perform other standard operations to the normal C/C++ variable as desired.

Note the use of pointer d* and non-pointer d1 is just to show that both ways are possible.


Using with Fortran
====================
Brook+, being like an extension to C/C++, is better called from C/C++ functions. But, provided that C/C++ wrappers are built for the Brook+ functions and then packaged into a DLL library, then any other language, eg Fortran, can use the GPU by calling on the C/C++ wrappers in the DLL.
Brook+ functions <--- C/C+ wrappers <--- Windows DLL / Unix shared objects <--- Fortran

Thursday, August 02, 2007

Notes Fortran

NotesFortran

Contents
=========

FORTRAN95 features
!DEC$ - compiler Directives
Keyword - SEQUENCE
Keyword - PURE
Initialize Data with slash /
Subroutines passed as argument of another subroutine
Allocatable Pointers / Array of Pointers
Using pointers in procedures
Dynamic Memory Allocation with Pointers
Declaration Statements for Arrays
Derived Data Types - with pointer components and being used as dummy pointer arguments
OpenMP
Static, Stack, Heap
Converting Integer to Character OR Writing to variable
Handling Character Strings between C# and Fortran
OpenMP programming warnings
Fortran Editors
String Manipulation



FORTRAN95 features
===================
List of F95 features implemented in Intel Fortran

1. FORALL -
2. PURE - for safety - ensure only INTENT(OUT,INOUT) arguments are changed.
3. ELEMENTAL - a type of PURE routine. Allow operation of arrays on element level.
4. CPU_TIME -
5. NULL intrinsic function - allow allocatable arrays to be pointed to this.


!DEC$ - compiler Directives
============================
Compaq (Digital) Visual Fortran (http://www.canaimasoft.com/f90vb/onlinemanuals/usermanual/TH_60.htm)


Compaq's Visual Fortran compiler (CVF) offers a great deal of flexibility to
create DLLs callable from C and/or Visual Basic. The compiler has a good set of
options to modify name-mangling, calling conventions and the method used to pass
arguments. Most of these options can be indicated through the use of compiler
directives (or pragmas) embedded in the source code. CVF compiler directives are
defined as comment lines, starting with DEC$.


To indicate that a subroutine should conform to the standard calling convention,
you add the DEC$ATTRIBUTE STDCALL compiler directive to the declaration of the
subroutine. For example:


subroutine MySub(argument1, argument2)  
   !DEC$ATTRIBUTES STDCALL:: GenDNASequence

In Compaq Visual Fortran, adding the STDCALL attribute also changes the default
method used to pass arguments, so you need to tell the compiler that arguments to
subroutine MySub are passed by reference. You can do this with the
DEC$ATTRIBUTE REFERENCE directive:

subroutine MySub(argument1, argument2)
!DEC$ATTRIBUTES STDCALL:: MySub
!DEC$ATTRIBUTES REFERENCE:: argument1, argument2


Also, when a procedure is declared with the standard calling convention (STDCALL),
Compaq Visual Fortran mangles its name. The name-mangling performed by CVF converts
the name of the procedure to all uppercase, adds an underscore as a prefix to the name,
and appends an at symbol (@) followed by the size of the stack (in bytes) at the end of
the name. The size of the stack is equal to 4 times the number of arguments in the
subroutine. MySub has 2 arguments, so the mangled name will be:


_MYSUB@8


Having to use this name to call our DLL procedure would be awful in most languages,
and illegal in Visual Basic. You can use another DEC$ATTRIBUTES compiler directive to
indicate an alias for the mangled name of the exported subroutine:



!DEC$ATTRIBUTES ALIAS: 'MySub'::MySub



The first argument is the alias (i.e. the name by which the subroutine would be available
to external programs using the DLL), the second argument is the Fortran name of the subroutine.



Finally, to indicate that the subroutine must be exported to the DLL as a public procedure,
you add the DEC$ATTRIBUTES DLLEXPORT compiler directive to the body of the subroutine:

!DEC$ATTRIBUTES DLLEXPORT:: MySub

So the full declaration of MySub would look like this:

subroutine MySub(argument1, argument2)
!DEC$ATTRIBUTES STDCALL:: MySub
!DEC$ATTRIBUTES DLLEXPORT:: MySub
!DEC$ATTRIBUTES ALIAS: 'MySub'::MySub
!DEC$ATTRIBUTES REFERENCE:: argument1, argument2
…


To create a DLL, you compile and link using the /dll switch. The following command-line
would compile and link MySub.f90 (containing subroutine MySub) into MySub.dll:


f90 MySub.f90 /dll /out:MySub.dll


Keyword - SEQUENCE
===================
SEQUENCE cause the components of the derived type to be stored in the same sequence they are
listed in the type definition. If SEQUENCE is specified, all derived types specified in component
definitions must be sequence types.


Keyword - PURE
===================
Pure Procedures
A pure procedure is a user-defined procedure that is specified by using the prefix PURE (or
ELEMENTAL) in a FUNCTION or SUBROUTINE statement. Pure procedures are a feature of
Fortran 95.
A pure procedure has no side effects. It has no effect on the state of the program, except for the
following:    
• For functions: It returns a value.
• For subroutines: It modifies INTENT(OUT) and INTENT(INOUT) parameters.
The following intrinsic procedures are implicitly pure:    
• All intrinsic functions  
• The elemental intrinsic subroutine MVBITS
• The intrinsic subroutine MOVE_ALLOC
A statement function is pure only if all functions that it references are pure.



Initialize Data with slash /
==============================
Variables are not auto-initialized. To initialize a variable as you declare it,
put the initial value between two slashes. This kind of initialization is done once
when the unit is first loaded, and hence, it should not be used in a subprogram that
gets invoked repeatedly (use an assignment instead). For symbolic constants,
initialization is achieved via the parameter statement.

    character title*20 / 'York' /
    integer*2 count / 0 /
    real*4 amount / 1.0 /
    !---------------------------------------------------
    ! Note that title will have 20 characters even though
    ! we stored only 4 (they will be padded by blanks).


Subroutines passed as argument of another subroutine
======================================================
subroutine subA( subB )

The subroutine named subB needs to be declared EXTERNAL and hence not be in a module.
If a function needs to be in the same module then can be done as follows:

module test
  contains
     subroutine subA(subB)
external subB
     .....
     subroutine subC()
end module
subroutine subB()
   use test
   call subC()
end subroutine


Allocatable Pointers / Array of Pointers
=========================================
REF: Fortran 90/95 for Scientists and Engineers, Stephen J Chapman
1. It is illegal to have array pointers of native data types in Fortran:
   REAL, DIMENSION(:), POINTER :: PTR
   - the dimension attribute refers to the pointer's target, not of the pointer.
   - the dimension must be deffered shape and the size is the size of the target, not
     the pointer.
2. It is legal to have array pointers by using derived data types:
  TYPE :: ptr
     REAL, DIMENSION(:), POINTER :: P
  END TYPE
  TYPE(ptr), DIMENSION(3) :: P1

REF: Intel Fortran Language Reference
1. In contrast to allocatable arrays, a pointer can be allocated a new target even if it is currently associated with target. The previous association is broken and the pointer is then associated with the new target.
2. If the previous target was created by allocation, it becomes inaccessible unless it can still be referred to by other pointers that are currently associated with it.


Using pointers in procedures
===============================
Pointers may be used as dummy arguments in procedures and may be passed as actual arguments to procedures.
1. If a procedure has dummy arguments with either POINTER or TARGET attributes, then the procedure must have an explicit interface.
2. if a dummy argument is a pointer, then the actual argument passed to the procedure must be a pointer of the same type, kind and rank.
3. a pointer dummy argument must not have the intent attribute
4. a dummy argument cannot appear in an ELEMENTAL procedure in Fortran95.


Dynamic Memory Allocation with Pointers
========================================
REF: Fortran 90/95 for Scientists and Engineers, Stephen J Chapman

REAL, DIMENSION(:), POINTER :: ptr1
ALLOCATE (ptr1(1:10))

This statement creates an unnamed data object of the specified size and the pointer's type and sets the pointer to point to the object. Because the new data object is unnamed, it can only be accessed by using the pointer. After the statement is executed, the association status of the pointer becomes associated. If the pointer was associated with another data object before the ALLOCATE statement is executed, then that association is lost.

The data object created by using the pointer ALLOCATE statement is unnamed and so can only be accessed by the pointer. If all pointers to that memory are either nullified or reassociated, with other targets, then the data object is no longer accessible by the program. The object is still present in memory, but it is no longer possible to use it => MEMORY LEAK.

If a piece of allicated memory is deallocated, then all pointers to that memory should be nullified or reassigned. One of them is automatically nullified by the DEALLOCATE statement, and any others should be nullified in NULLIFY statements.

Declaration Statements for Arrays
===================================
SUBROUTINE SUB(N,C,D,Z)
   REAL, DIMENSION(N,15) :: IARRY        !explicit shape array
   REAL, C(:), D(0:)                     !assumed shape array
   REAL, POINTER :: B(:,:)               !deferred shape array pointer
   REAL, ALLOCATABLE, DIMENSION(:) :: K  !deferred shape allocatable array
   REAL :: Z(N,*)                        !assumed size array

Automatic Arrays - local array in a function, whose size is one of the arguments
Adjustable Arrays - array which is an argument, whose size is also one of the arguments.

To use arrays efficiently, see Optimizing Applications -> Programming Guidelines -> Using Arrays Efficiently.
When passing arrays as arguments, either the starting (base) address of the array or the address of an
array descriptor is passed:
When using explicit-shape (or assumed-size) arrays to receive an array,
the starting address of the array is passed.

When using deferred-shape or assumed-shape arrays to receive an array,
the address of the array descriptor is passed (the compiler creates the array descriptor).

Automatic vs Save
Automatic variables can reduce memory use because only the variables currently being used are allocated to memory.

By default, the compiler allocates local variables of non-recursive subprograms, except for
allocatable arrays, in the static storage area. The compiler may choose to allocate a variable in
temporary (stack or register) storage if it notices that the variable is always defined before
use. Appropriate use of the SAVE attribute can prevent compiler warnings if a variable is used
before it is defined.




Derived Data Types - with pointer components and being used as dummy pointer arguments
========================================================================================
Given a derived data type with pointer arguments:
   type varDDT
      type(BigDDT), pointer :: rComp
   end type
Given that it is used as a dummy pointer argument
   subroutine foo(aDDT)
      type(varDDT), pointer :: aDDT
 ....
    end subroutine

Then, in another function which uses "foo", the DDT can be used as:
program
  type(varDDT), pointer :: aDDT_p
  type(varDDT) :: vDDT

  call testFoo(vDDT)   ! vDDT contains the data found in testFoo
  nullify(aDDT_p)
    end program

    subroutine testFoo(vDDT_local)
  type(varDDT) :: vDDT_local
       aDDT_p => vDDT_local

  call foo(aDDT_p)
       ! vDDT_local will be able to pass to outside routine safely.
    end subroutine


Note that this method is not required if the DDT concerned does not contain components
which are also DDT.


OpenMP
=======

Prerequisite:
Before inserting any OpenMP parallel directives, verify that your code is safe for parallel execution by doing the following:

Place local variables on the stack. This is the default behavior of the Intel Fortran Compiler when -openmp is used.

Use -automatic (or -auto_scalar) to make the locals automatic. This is the default behavior of the Intel Fortran Compiler when -openmp is used. Avoid using the -save option, which inhibits stack allocation of local variables. By default, automatic local scalar variables become shared across threads, so you may need to add synchronization code to ensure proper access by threads.

Static, Stack, Heap
====================
Heap area stores dynamic arrays
Static storage area store variables that are available for the life time of the program.
In C, local variables are stored in the stack.

In Fortran,
"By default, the compiler allocates local scalar variables on the stack. Other, non-allocatable variables of non-recursive subprograms are allocated in static storage by default. This default can be changed through compiler options. Appropriate use of the SAVE attribute may be required if your program assumes that local variables retain their definition across subprogram calls."

For openMP, local variables in Fortran will become automatic and thus stored in the stack. Note that
OpenMP threading model is based on threads and each has its own stack.



"static" as a descriptive term refers to the lifetime of C++ memory or storage locations. There are several types of storage:
        - static
        - dynamic (heap)
        - auto (stack)
A typical storage layout scheme will have the following arrangement, from lowest to highest virtual memory address:
        text (program code)
        static (initialized and uninitialized data)
        heap
        (large virtual address space gap)
        stack
with the heap and stack growing toward each other. The C++ draft standard does not mandate this arrangement, and this example is only an illustration of one way of doing it.


heap-array in Intel Fortran - This option puts automatic arrays and arrays created for temporary computations on the heap instead of the stack.
on Windows:
    /heap-arrays-        = no heap arrays (default)
    /heap-arrays[:size]  = where arrays of size (in kb) or larger are put on the heap.
on Linux:
    -no-heap-arrays        = no heap arrays (default)
    -heap-arrays [size]    = where arrays of size (in kb) or larger are put on the heap.
Example:
    In Fortran, an automatic array gets it size from a run-time expression. For example:
RECURSIVE SUBROUTINE F( N )
INTEGER :: N
REAL :: X ( N )     ! an automatic array
REAL :: Y ( 1000 )  ! an explicit-shape local array on the stack Array X in the example above
                    ! is affected by the heap-array option. Array Y is not.





Converting Integer to Character OR Writing to variable
========================================================
Example:
   CHARACTER(LEN=15), ALLOCATABLE  :: cPctile(:)
   ALLOCATE( cPctile(nPctile) )
   write(cPctile(jj), '(F15.7)') pctile(jj)

The number in pctile(jj) is being written into a CHARACTER variable called cPctile. Note the character has length of
15 which is enough to write 15 characters specified by the Format F15.7.



Handling Character Strings between C# and FortranSubmit New Article
Last Modified On :   December 6, 2009 6:09 PM PST
Rate Please login to rate! Current Score: 0 out of 0 usersPlease login to rate! Current Score: 0 out of 0 usersPlease login to rate! Current Score: 0 out of 0 usersPlease login to rate! Current Score: 0 out of 0 usersPlease login to rate! Current Score: 0 out of 0 users






Handling Character Strings between C# and Fortran
=================================================================

Passing Character Strings as In parameters from C# to Fortran

C# provides built-in reference type "string" representing a string of Unicode characters. It is an alias for String in the .NET Framework. When passing "string" type by value from C# to Fortran function or subroutine platform invoke service copies string parameters, converting them from the .NET Framework format (Unicode) to the unmanaged format (ANSI), if needed. The unmanaged format is null-terminated so the C# method prototype of Fortran function or subroutine must account for the length argument passed along with the string address.


Passing Strings as In/Out parameters from C# to Fortran

Managed strings are immutable, platform invoke does not copy them back from unmanaged memory to managed memory when the function returns. If the Fortran function or subroutine wants In/Out parameters you need use StringBuilder Class in C#. StringBuilder Class represents a string-like object whose value is a mutable sequence of characters.

Fortran subroutine
subroutine FPassStringSub (Int_Arg, Str_In, Str_Out)
!DEC$ ATTRIBUTES DLLEXPORT :: FPassStringSub
integer, intent(in) :: Int_Arg
character*(*), intent(in) :: Str_In
character*(*), intent(out) :: Str_Out
end subroutine



C# method prototype
        [DllImport("FDLL.dll", CharSet = CharSet.Ansi, CallingConvention = CallingConvention.Cdecl)]
        public static extern void FPASSSTRINGSUB(ref int Int_Arg, string Str_In, StringBuilder Str_Out, int Str_In_Len, int Str_Out_Len);
STR_IN_LEN, int STR_OUT_LEN);



Returning Character Data Types from Fortran to C#

Similar to how C language handles Fortran function returning character data type the corresponding C# method prototype must add two additional arguments to the beginning of the argument list:
-  The first argument is a StringBuilder object where the called function should store the result.
-  The second argument is an int value showing the maximum number of characters that must be returned, padded with blank spaces if necessary.


Fortran function

function FRetString(Int_Arg)
!DEC$ ATTRIBUTES DLLEXPORT :: FRetString
character(*) :: FRetString
integer, intent(in) :: Int_Arg
end function


C# method prototype
        [DllImport("FDLL.dll", CharSet = CharSet.Ansi, CallingConvention = CallingConvention.Cdecl)]
        public static extern void FRETSTRING(StringBuilder Str_Result, int Res_Len, ref int Int_Arg);



How to avoid Race Conditions in OpenMP programming
===================================================
This section aims to cover a few points to watch for when doing parallel programming using OpenMP in Fortran.
The tips here could also apply to different languages with OpenMP such as C/C++.
1. First danger sign of parallel programming bug is inconsistent behaviour.
2. When a program stalls or hangs indefinitely sometimes, but other times run to completion.
3. When a program crashes sometimes, but other times run to completion.
4. When OpenMP is applied to a higher level function that calls many small low level functions, check  whether module variables have SAVE attribute.
5. If module variables have save attributes, then check if those variables can be declared PRIVATE too so that no outside functions can use them.
6. When the SAVE module variables are being altered in a function in the module, check to see if the function can be encapsulated in a OMP parallel region.
7. When checking for race conditions using tools like Intel Thread Checker, ensure that the number of threads specified is not greater than the number of processors / cores.


Fortran Editors
=================

The following list are GUI based Fortran editors development environment.

Force
http://force.lepsch.com/

Plato - for Silverforst FTN95
http://www.silverfrost.com/16/ftn95/plato.aspx

Photran - on Eclipse IDE
http://www.eclipse.org/photran/

Geany - GTK based
http://www.geany.org/Main/AllFiletypes


There are many other general purpose editors that can edit Fortran including: Vi, Emacs, EditPlus.

String Manipulation
====================
There is a misconception that the string manipulating and handling capabilities are limited. However, it may not be as limited as initially thought. Here are a few useful string handling features from poplular languages that also exist in Fortran.

Comparing Substring:
   INDEX( string, substring) -> returns integer representing the first position of occurance of the substring in string.
   Converting Integer to String
write(charWord, '(I5)') ii-1
where ii is an integer which is written to the character variable 'charWord'.