今日の気になる曲。

• 
  9:30 切れの直前、（なよなよ系）男性ボーカル、特徴的なフレーズのリピート。
• 
  付記： Iron & Wine の White Tooth Man という曲らしい。 （youtube にあったので、確認済み。）

ryuon

• 
  nonlinear solver の組み込み。
  • 
    ふと思い立って、以前からやっておかなければと思っていたことに再びアタックする。
    • 
      Newton-GMRES に関して以前調べたとき （11/23/2007）に見つけた NITSOL という library を使ってみようという話。
    • 
      NITSOL の reference は以下のとおり：
      M. Pernice and H.F. Walker, SIAM J.Sci.Comput. (1998) 19, pp.302-318:
      "NITSOL: A Newton Iterative Solver for Nonlinear Systems"
      (http://www.siam.org/journals/sisc/19-1/30384.html)
  • 
    FORTRAN のコードを C から使う、という毎度ある事例で手こずるのだが、何とか出来た。
    • 
      以前もはまって、解決したような気がしていたが、 その時の記述などが見当たらないので、また１から調べ直した。
    • 
      ポイントを「FORTRAN のライブラリを C から使う事例の１つとして」にまとめておく。
  • 
    libstokes に登録しておく。
    • 
      まだ top level （例えば stokes3 の level）から 使えるようにはなっていない。

FORTRAN のライブラリを C から使う事例の１つとして。

• 
  FORTRAN の common block の C からの使い方：
  • 
    FORTRAN の common block

    
    
          common /nitprint/ iplvl, ipunit
    
       

    へのアクセスは、構造体

    
    
    extern struct {
    
      int iplvl;
    
      int ipunit;
    
    } nitprint_;
    
       

    を介して行える。
  • 
    使い方の例 in C：

    
    
      nitprint_.iplvl = 1;
    
      nitprint_.ipunit = 6;
    
       

• 
  引き数に関数がある FORTRAN のルーチンの呼び方、および FORTRAN で使われる関数の C での実装：
  • 
    難しいと思い込んでいたが、普通に書いたら普通に出来た。
  • 
    FORTRAN の subroutine

    
    
          subroutine nitsol(n, x, f, jacv, ftol, stptol, 
    
         $     input, info, rwork, rpar, ipar, iterm, dinpr, dnorm)
          implicit none  
          integer n, input(10), info(6), ipar(*), iterm 
    
          double precision x(n), ftol, stptol, rwork(*), rpar(*) 
    
          double precision dinpr, dnorm
    
          external f, jacv, dinpr, dnorm
    
       

    というものがあった場合、 つまり void 型の関数 f と jacb と double を返す関数 dinpr と dnorm があって、それぞれ FORTRAN では

    
    
          subroutine f(n, xcur, fcur, rpar, ipar, itrmf)
    
          implicit none
    
          integer itrmf, n, ipar(*)
    
          double precision xcur(n), fcur(n), rpar(*)
          subroutine jacv(n,xcur,fcur,ijob,v,z,rpar,ipar,itrmjv)
    
          implicit none
    
          integer n, ijob, ipar(*), itrmjv
    
          double precision xcur(n), fcur(n), v(n), z(n), rpar(*)
          subroutine dinpr( n, x, sx, y, sy )
    
          integer n, sx, sy
    
          double precision x, y
          subroutine dnorm( n, x, sx )
    
          integer n, sx
    
          double precision x
    
       

    と定義されている場合を考える。
  • 
    （例１） C から nitsol だけを呼びたい場合 （f と jacv は FORTRAN のものを使う）
    • 
      FORTRAN の関数の C での function prototype declarations:

      
      
      void nitsol_ (int *n, double *x,
      
      	      void (*f)(int *n, double *xcur, double *fcur,
      
      			double *rpar, int *ipar, int *itrmf),
      
      	      void (*jacv)(int *n, double *xcur, double *fcur,
      
      			   int *ijob, double *v, double *z,
      
      			   double *rpar, int *ipar, int *itrmjv),
      
      	      double *ftol, double *stptol,
      
      	      int *input, int *info, double *rwork, double *rpar,
      
      	      int *ipar, int *iterm,
      
      	      double (*dinpr)(int* N, 
      
      			      double* X, int* incX, 
      
      			      double* Y, int* incY),
      
      	      double (*dnorm)(int* N, double* X, int* incX));
      void fbratu_ (int *n, double *xcur, double *fcur,
      
      	      double *rpar, int *ipar, int *itrmf);
      void jacvbratu_ (int *n, double *xcur, double *fcur,
      
      		 int *ijob, double *v, double *z,
      
      		 double *rpar, int *ipar, int *itrmjv);
      // BLAS functions
      
      double
      
      ddot_(int* N, 
      
            double* X, int* incX, 
      
            double* Y, int* incY);
      
      double
      
      dnrm2_(int* N, 
      
             double* X, int* incX);
      
          

    • 
      C からの呼び出し：

      
      
        int n = MAXN;
      
        double *x = (double *)malloc (sizeof (double) * MAXN);
      
        double *rwork = (double *)malloc (sizeof (double) * LRWORK);
      
        double *rpar = (double *)malloc (sizeof (double) * LRPAR);
      
        double ftol = 1.0e-6;
      
        double stptol = 1.0e-6;
      
        int input[10];
      
        int info[6];
      
        int ipar[2];
      
        int iterm;
        nitsol_(&n, x, f_, jacv_, &ftol, &stptol, 
      
                input, info, rwork, rpar, ipar, &iterm, ddot_, dnrm2_);
      
          

  • 
    （例２） C で f を定義したい場合

    
    
    void f_c (int *n, double *xcur, double *fcur,
    
    	  double *rpar, int *ipar, int *itrmf);
    
         

    という関数を書いて、以下のように普通に呼び出す。

    
    
      nitsol_(&n, x, f_c, jacv_, &ftol, &stptol, 
    
              input, info, rwork, rpar, ipar, &iterm, ddot_, dnrm2_);
    
       

• 
  すべて GCC と GNU Fortran compiler の話。
  • 
    他の環境は知らない。といっても、実際上何があるのかな？ Intel の compiler くらいかな、思いつくのは。
• 
  references:
  • 
    6/4/2008: [ryuon] NITSOL の組み込み。
  • 
    5/8/2003: 行列に関する(個人的な)混乱。

帰宅したらオランダから手紙が届いていた。

• 
  あちらが送ったという知らせを受けたのが 5/15 だったから、 結構かかったことになるな。やっぱり遠い（のか、カナダの郵便屋さんの問題なのか……）。