fortran tips と自己流書き方ガイド

今のfortranは良くなった。コンパイラバグが激変したし、aptでmklをインストールできblas,lapack,fftwなどがおよそ標準装備となった。そもそも数学関数が高速になりいちいちblasを呼ばなくていいようになったし、ベクトル表記や関数などがだいぶと整備されてきた。もちろんおかしなところも多々あるとも言えるが古いコードが動くのも良い
（python2-->python3ではひどい目にあった。いまだにsystem関係がメタメタでsubprocessの仕様はひどくていずれ置き換えられるだろう。os.systemしか信用してない）。

良いコードを書くための自己流ガイドライン

（あくまで柔軟対応が基本。経験は少ししか役立たない）。

• 理想はコメント行の不要なコード あるいは論文の数式との対応だけが書いてあるだけで理解できるようなもの。コメントを読まないとわからないようなプログラムは問題がある。暗号のような符牒を工夫しないほうがいい---計算アルゴリズムの本質と無関係な論理構造は極力持ち込まない。「便利な関数化」をしたつもりが保守性の悪い「ブラックボックス\inブラックボックス\inブラックボックス...」にならないように。
• 情報の二重記述をしないようにする. たとえば本体と呼び出し側での整合性をプログラマに要求するような書き方は極力避けたい。「二重記述がダブルチェックになる」というようなことはほとんどない。極力interface文は書かないー自動決定されるような部分を人間が書くべきではない。moduleを活用するのはインターフェイス文を自動で作ってくれるという意味においても有効。たとえばpure関数,optional変数やsize関数の有効化、名前結合、などをさせるには、関数やサブルーチンをmoduleに入れてやりuseで用いる必要がある。
• module変数は（ごく少数で明示的な例外を除いては）protectすること。デフォルトをprivateとしpublicを書く。useするときは必ずuse,onlyで明示する。
• 変数の居場所（どのmoduleに属するか？）を決めること。データの唯一性を守ること。常にその変数の属するmoduleをuseして読み出す。たとえば原子数を指す変数は唯一でないといけない。そのコピーを作ってはならない。（極端で変な例だけど：データの入った行列の本体aがmodule Aに、そのサイズを指定する変数nがmodule Bに入っているのは問題がない。問題なのはそのサイズ変数nをmodule Aとmodule Bに重複して入れてしまうこと。）。また、たとえばある場所で、突然電子密度なり、原子基底のインデックスなりが必要になったりする、このときモジュールであれば、use onlyのみで済む。サイドエフェクトなく引用できる。変数を増やすなどのコードの改変も本質的に無関係なところをいじる必要がない。「データを読むのはどこからでも自由に読める（onlyで明記できる.変なデータ経路を作らずに直接参照で書くべし！）、生成するのは一箇所でのみ」を守る。
• pytorchなどでは、tensorは必ず次元・サイズの概念と対になっている。こういうやり方もありだと思う。しかし、fortranでは、次元・サイズの自動結合のチェックはpytorchほどうまくいかなくて面倒になりうる。どうしてもgemmなどを使うことも必要なので。
• 要するにモジュールを唯一性クラスオブジェクト（シングルトン）として利用することになる。複数作れないので不便だと思うかもしれないけれどルールが単純になるご利益もある；唯一性があるので名付ける必要がなくプログラム起動時から存在するインスタンスとなる。selfがいらないし、唯一性のため混乱を生じずどこからでもダイレクトに読める。具体的には
  • まず「ファイル読み込みはモジュールでおこなって(allocateもおこなって）そのデータはprotectしておいて他から読み出すときは必ずuse onlyする」とするのがよい。
  • 初期化はbootstrap方式で連鎖的にモジュール初期化ルーチン（引数ほとんどなし）を呼んで、複数のモジュールに初期化データを貯めていく。これらの変数は計算が終了するまで変更されないワンスライトなデータセットである。
  • サブルーチンはモジュール内に書く。サブルーチンが必要とする入力データはuse onlyでどこかのモジュールから読みだす。自分自身で言えば、昔のサブルーチンbndfpは引数が50個とかあったんだけど、計算指示のための引数3個ぐらいに削減しました。
  • モジュール名はm_foobarとする。
  • 計算の流れをコントロールする指示的な変数のみサブルーチン入力引数で渡す。ただ、フェルミエネルギーとか、キーになる物理量を指示引数として取り扱うのはありえるー計算フローが読みやすくなるメリットがありえる。
  • 出力値はモジュール変数で戻す（サブルーチン内でallocateする）のを基本とする。サブルーチン引数を出力値としない。ただし、そうするとuseする変数が出力値となるから「出力値が関数コールから離れたところに書かれている」ということがおきる。ただcyclic searchのできるエディタであれば対応できる。メインルーチンにおいてサブルーチンAのコールでモジュール内変数aが変更されそれを次のサブルーチンB内で使うということも発生しうる点、すなわちデータ生成の手順が見えにくくなりうる点があり注意は必要ー明瞭になるようにしておく必要がある。
  • 要するにサブルーチンの仕事は「複数の必要なモジュール変数を読み取って、なにか計算してそのサブルーチンを含むモジュールのモジュール変数をallocateしてデータを書き出す」という仕事となる。そのデータは唯一性を持つものとして扱いそのデータが必要なときにはダイレクトにそれを読むこと。サブルーチンの引数で渡さない。
  • 小さなサブルーチン（基本的にI/Oの明瞭な小さい関数や数学関数など）は独立関数とする。あるいはモジュールに格納しておくー*.modというインターフェイスを自動生成してくれるので最低限の型チェックをコンパイラに任せることができる。モジュール変数はなくて良い。最新のgfortranだと-fallow-argument-mismatch
    をつけないと古いコードは動かなかったりするーこれを考えてもモジュールに格納する方式にしていくのがよい。たとえばmodule m_utilというモジュールを作ってcontainsで多数のサブルーチンを格納する。

ぐらいに注意すればなんとかなる。

• そもそもfortranの限度をわきまえるべき。fortranで書けるのはfile to fileのfeedfowardなプログラム（ほとんどすべての変数はワンスライトなもの）ぐらいに心得ておく。グルーランゲージで外側から小さなfortranプログラムを連鎖的に実行させる。iterationしたいなら外側で回す。現状、自身の電子状態計算のコードではそうなっていないが将来的には（もうすこしfortran+pythonが進化すれば）それが良いと思っている。
• 将来予測はあまり役に立たないが少しだけ考えている.ぼくの計算では複数のfortranプログラムをpythonスクリプトos.systemで順に呼び出しイテレーションサイクルで計算させている。fortranプログラムのメインルーチンは、形式的にはサブルーチン（引数なし）でbind(C)してあり、ライブラリ化している。pythonのmpi4pyからcommをfortranに渡す仕組みをテストしてみたが現状では面倒なことが多いので,pythonをメインルーチンにする方法をデフォルトにはできてない。python＋fortranならよいが、python+C+fortranとしたいとは思ってない。
• block文はとても有効。単なるラベルをつけてのブロック区切りとしても使える。延々と長いコメントに似たラベルでも良い。ヘタにサブルーチン化するまえにブロック化すればよい。if文なども含めてラベルをしっかり書く。局所変数やuseを使えるしallocateをかなり減らすことができる。emacsなど自動でendblockにラベル補完してもらえるエディターを使う。
• インプット変数の多いサブルーチンはmodule化する。ほとんどのインプット変数は「環境変数」的なもので、初期化時点でどこかのモジュールに属させれば良いものであったりする。
• moduleに入れたサブルーチンをuseして使う場合にはコンパイラが変数の型整合性をチェックしてくれる。なのでサブルーチンは（モジュール変数がない場合でも）モジュールに入れてuseして使ったほうがいい。並列コンパイルmake -jでやればまあまあ速度は何とかなる。
• logical,integer,real(8),complex(8),character(バイト数）を基本にする。kindなどと書かない、と思っていたが混合精度をやり始めるとどうしても使わざるを得ない点はある。まあごく部分的に導入すればいいのだけど。
• 文字列操作はそれなりにできる。が、無理がたたりそうなら,call system(foobar)でpythonでも呼ぶ。ファイルtoファイルでやり取りする。何千行かあったfortranコードによる文字列操作部分をcall system(foobar.py)に直したら100行程度になった（と同時に（問題にならない部分ではあったが）pythonの遅さを再認識した）。
  -ベクトル的表記に置いて（不整合が起こりえるだけでメモ的な意味しかない場合）、配列サイズは(pytorch方式のように）コメントとしてメモ書きするだけのほうが良かったりもする.コーディングは（なんとなく、ではなく）、pros and consを考えて行うのが理想（文章もそうだけど）。ただ最善手を一意的に決めるのは困難なときも多い。
• いまだにopen(newunit=...)を紹介していないfortran学習文献もある。マニュアルを見たほうがよいことも多い。新機能は実装できてないこともあり得るし実装されててもコケることもあるのでgfortranとifortの二種のコンパイラは通す（スタンダード実装であるがfortranの正式な取り決めにない拡張もありえる）。正直、英語圏と日本語圏では情報量に差があるなと感じることも多いー（できたら避けたいけれども）英語検索も必要。調べているとひょっとすると「日本語だと大多数の教育用文献でnewunitを紹介していない。」のかもしれないー初心者はマニュアル見ないだろうし不安になってきました。今はbingやgithubのcopilotに聞く、という手もあるがある程度の理解がないとそもそも聞けないしウソかホントか判別できないというのもある。
• サブルーチン化、モジュール化、クラス化に本当に意味があるか？本質的に論理的なクリーンさを与えることができているのか？を考える。脳の無駄遣いになるような新規概念の導入はやめよう。特に構造体＝＞クラス化は柔軟性に乏しく危険な気がする。サブルーチン化することできれいになってるようでなっていないことも多い。本質的に、一度しか呼ばないサブルーチン、というのは意味がない。モジュールあるいはblockで良い。
• CMakeを使うべき。結局はCMakeFiles.txtを書くのが良い。github copilot chatはまあまあよく教えてくれるのでCMakeなんか勉強したくない、という人でもCMakeFiles.txtを書けるようにはなってきている。
• どうしても古いコードとの共存になる。とにかくモジュールに入れていくのがいいかもしれない。指示変数（サブルーチンの冒頭で渡す）、モジュール変数（モジュールの出力）、use onlyによる読み込み、に分ける。まずは初期設定的なサブルーチンをモジュール化し、その変数をuse onlyするようにする。
• 必ずエンバグするので、テストシステムを工夫しておかないといけない。テストシステムとgitを使わないと開発が維持できない。昼飯とかのスキマ時間にこまめに起動する。
• コンパイラバグがないわけではない。-O0で通るならコンパイラバグの可能性が高い。-O0でもあかんのならたぶんコンパイラバグではなかろう。どのソースを-O0にするかを見つけて、conditional compileの指示をmakefileに書き込まないといけない。
• バグ取りにはオプション(-gとか-tracebackとか）をつけて、問題ありそうなコードをコンパイルし直して実行する。まずはシングルでも落ちるかどうかを試す。デバッガは面倒だし役に立たない（というかぼくが使えてないだけか）。数値的エラーだとチェックライトをいれるしかないし、基本的にはバグを追い詰めていくような方法しかない。バグは複数ありえるので分割的に追い詰める。ホームズにでもなった気分で追い詰める。論理力が試されている。まじめにやればいい、というのものでない。最短時間で最大効率を上げるデバッグ法をその都度その都度発明しないといけない。エンバグの場合、gitを活用する必要もある。
• 混合プログラミング。どうしようもない場合以外、いたずらな複合プログラミングは避ける。へたにbind(C)などとするのは面倒になりうる。call systemでファイルtoファイルのデータIOでやれるならそれがよいかもしれない。pythonを親にすえるのなら、ごく限定的なインプット引数のみ（あるいは引数なし）のfortran関数のコールがいいかもしれない。独立したプログラムを連鎖的にコールしていくのと似てはいるが、モジュール内部の変数を保った状態で、連鎖的にコールしていけるアドバンテージはある。
• 整数化するには、nint,ceiling,floorを使い分ける。負数に対してintは定義が面倒。
• 以下にいくらか示すがfortran intrinsic functionsには結構便利なものがある。数値計算の多くはindex計算だったりする。pytorchほどの自由はきかないが極力ベクトル演算とし、do文、深い階層化は極力避ける。
• 大きなdo loopにはラベルだけでなく数字もつけたほうが検索しやすいし視認性もあがるしコミュニケーションもとりやすい。
• subroutineのなかにcontainsでsubroutineを入れ込める機能がある。これは変数の受け渡しやスコープをメチャクチャにできるヤバイ機能。一見便利ですけどやめときましょう（ぼくは引数ゼロの書き出しだけのsubrouitne内subroutineを使っているけれどそのうち廃止したい）。
• fortranへの標準入力read(,)とかread(5,*)は撲滅しました。getargでいい。
• intent(inout)な関数は徐々に撲滅する方向で
• pointer(変数のエイリアス）は活用できる。doの中などで部分配列に対してポインターをあててやるのも良い。
• HDF5とかXMLとかJSONとかの階層化データフォーマットをfortranで読み書きするのはやりすぎな気がする。そもそもそういうのはだいたいすたれて難儀する。階層化させると柔軟性が乏しくなるしハンドリングも不必要に面倒になる。fortranのunformattedで書けば良い。direct accessも使える。人間が目視確認するところ、gnuplotへ出力するところなどだけformatをつける。mpiioは利用している。
• 一気にすべてを直すわけにも行かないので進化のさせ方が一番難しいところ。優先順位を考えないといけない。ぼくの場合は、新たな書き直しを半分試験的に埋め込んでいき、しばらく経って（計算に利用して）、問題なければ全面的にその書き直しを採用していく、というようなことをしている。
• coarrayは使ってない。歴史的にMPIで動かしてきたのであえてcoarrayに移行しないと思う。MPIとOpenACCを利用している。
• 自分が無知であることを認識する（自戒です）。過去のプログラム履歴を引きずっているのが普通なのですんなりとはいかない。教条主義もよくない。新たに学ばないといけない。それなりの頑固さも必要.それなりの柔軟さも必要.

sample collections

findloc

dim=1書かないと怒られる(gfortran),maxlocなどもある。
多くの部分でindex countを行っているのでこういうのは結構役立つ。論理式書けるので一般的なことができる。

ik3= findloc(kk3,value=kkk3(3),dim=1)

nn2 = findloc([(rsmh2(i,j)>0d0,i=1,nnx)],value=.true.,dim=1,back=.true.)

ispec(j)= findloc([(trim(alabl)==trim(slabl(i)),i=1,nspec)],dim=1,value=.true.)

iout= findloc([(a(i)>0,i=1,iin)],dim=1,value=.true.,back=.true.)-1

昔はこういうのなかったからクソループ書いてました。先生のコードにはクソループいっぱい入ってると思いますよ。もちろんぼくのコードにも。

reshape

zmelc = reshape(dconjg(zmel),shape=shape(zmelc),order=[3,2,1])

zsec= zsec+ matmul(reshape(zmelcww,[ntqxx,nm]),reshape(&
      reshape(zmel,shape=[ns2-ns1+1,ngb,ntqxx],order=[2,1,3]), [nm,ntqxx]))

GPUでもreshapeなんかはサポートしているようだ。再配置（order）などもシステムに任せちゃったほうがいいだろうと推察している。（昔は、subroutineコールで 引数をa(1,1)とするのとa(:,:)とするので驚異的に速度が違ってビックリしたこともあった。いまはmoduleにいれてるし、use
で渡しているのでそういう点にあまり気を使っていない）。

count

nt0  = count(ekc<ef) 

noccx(0:3) = count( ek_(1:nband, 0:3)<efermi,dim=1)

merge

merge,where,count,findlocなどを用いてif文やdo loopはできるだけ書かないこと。

 itini= merge(max(ns1,nt0m+1),  ns1, mask= omg>=ef)

source

mmm='abc'
allocate(farg,source=trim(mmm))

concurrent,forall

シンタックスが違う。

do concurrent (i2 = 1:mnl(ic),i1 = 1:mnl(ic)) 
...
forall(itp=lbound(zsec,1):ubound(zsec,1)) zsec(itp,itp)=1d0)

associate

ポインタ宣言いらなくて値を入れれる。以下のように数式でも可能。
ただ、変数だとポインタ的に振る舞うようだ。
https://www.nag-j.co.jp/fortran/fortran2003/Fortran2003_8_21.html

       do  io1 = 1, norb; if(blks(io1)==0) cycle 
          associate(k1=>ktab(io1),nlm11 => ltab(io1)**2+1)
            sumh = sumh+sum([( &
                 sum(qhh(k1,k2,ilm1,nlm21:nlm22,:)*ppihhz(k1,k2,ilm1,nlm21:nlm22,:)) &
                 ,ilm1=nlm11,nlm11+blks(io1)-1) ])
            do  ilm1 = nlm11, nlm11+ blks(io1)-1
               if( nlm21<= ilm1 .and. ilm1<=nlm21+blks(io2)-1) then
                  sumt = sumt + sum(qhh(k1,k2,ilm1,ilm1,:)*tauhh(k1,k2,ll(ilm1),:))
                  sumq = sumq + sum(qhh(k1,k2,ilm1,ilm1,:)*sighh(k1,k2,ll(ilm1),:))
               endif
            enddo
          endassociate
       enddo

配列構成子

便利なので多用している。sumなどと組み合わせることも多い。

以下のように多次元でも使える場合がある(gfortran,ifort)。matmaは自作の行列掛け算の外部関数。

         call matma(zmelcww,[(transpose(zmel(:,:,itpp)),itpp=lbound(zmel,3),ubound(zmel,3))],zsec,&
              size(zsec,1), size(zmelcww,2)*size(zmelcww,3), size(zsec,2)) !  zsec accumulated

fortranの理屈から言えば使えそうではある。ただ、いまどきのfortranなら怒ってくれてもいいようにも思うー外部関数は徹底的に関知しません、ってことか。（ちょっと非合法っぽいので使用はやめました）。

構造体の利用はむずかしい。

変数の分類学に困り放棄した。

一つの箱にまとめると便利と思うかもしれないが、未来はわからないのでまとめ方がどんどん変わりえる。データ構造は変わり得る。トンカチだけ使いたいのにいつも道具箱を持ちあるかないといけない、今日は薬もついでに入れとくか、ということにもなる。とちゅうで店に寄ってドライバを買って箱に入れるということにもなる。他者（あしたの自分）は「あ、薬が道具箱に入っている！」ということになる。じゃあっていうので薬箱と道具箱を作ってそれを収納する箱を作るとなる…　

また、インデックスシステム（たとえば原子だと、インデックスl,m,nの取りうる値が複雑に依存している）とそれに立脚するデータを別のモジュール（タイプ）に入れたいが、構造体を利用するとなるとムダな混乱が生じる。

電子状態計算では「初期的な環境的変数をもとにbootstrap式に、かなり多様な計算を行っていく」ので、まったくもって不便になりうる。

構造体のインスタンスはサブルーチンの引数として授受することになる。そのようなことをすると不要に長いデータ授受の経路が発生する。その経路を追う、ということが必要になる。もしそのようにせずに構造体のインスタンスを特定のモジュールに入れておいて読み書きするのなら、そもそもそのモジュール変数としておくだけでよい。箱の中の箱に入れなくて良い。結局のところ、構造体には束ねて隠すという本質しかない。関数も入れてクラス化できるというわけだが、おそらくモジュール関数（use onlyでデータにダイレクトにアクセスできる）で間に合うしスッキリする。

pointer配列

ただ、原子ごとに3種の動経方向の電子密度があるような以下場合では、以下のs_rv1のような汎用構造体を利用することができる。構造体をdeallocateするとその中身もメモリリークなく消してくれるようだ（たぶん）。原子ごとに違うメッシュサイズのデータを一括してコントロールできる。他の構造体はモジュールに書き換えたが、このpointer配列にのみ残っている。

use m_struc_def
type(s_rv1) :: rhoat(3,natom)
...
do i=1,natom
  allocate(rhoat(1,i)%v(1:nr(i))) !nrは動経方向のメッシュサイズ
  allocate(rhoat(2,i)%v(1:nr(i)))
  allocate(rhoat(3,i)%v(1:nr(i)))
enddo
使うときは
associate(v1=>rhoat(1,i)%v)
associate(v2=>rhoat(2,i)%v)

...
!========================
タイプの宣言
module m_struc_def
  type s_rv1
     real(8),allocatable:: v(:)
  endtype s_rv1
...

intentは桁を揃えて書く。

変数宣言とは別にする。use やらimplicitを先に書かないといけない仕様は変だけどしかたない。

module m_rotwave
  public:: rotmto,rotmto2,rotipw,rotipw2,rotwvigg
contains
  !!------------------------------------------------------
  subroutine rotwvigg(igg,q,qtarget,ndimh,napw_in,nband,evec,evecout,ierr)
    intent(in)::      igg,q,qtarget,ndimh,napw_in,nband,evec
    intent(out)::                                            evecout,ierr
...

（としたいが、実際にはuseとimplicit noneの下へ持っていかないと動かない）

(いずれは出力変数はpythonの関数のように左辺に書けるようなるだろうけど。
関数とサブルーチンを区別しなくてもいいだろうし）。

できるだけサイズを渡さない

基本的にモジュール内でsubroutineを定義しているのなら、変数を渡せば、
変数のサイズ情報も受け渡すことができる。

なので以下が機能するには、gtv2_setrcdをモジュール内で定義し、use gtv2_setrcdしてから使う必要がある。

nrecsなどはモジュール変数。
...
  subroutine gtv2_setrcd(recrdin)
    character(len=*):: recrdin(:)
    nrecs = size(recrdin)
    reclnr= len(recrdin(1))   !write(6,*)nrecs,reclnr
    allocate(character(reclnr)::recrd(nrecs))
    recrd=recrdin
  end subroutine gtv2_setrcd

moduleの使い方

• 以下では初期化m_igv2xall_initでモジュール変数にデータをセットする。
  関数m_Igv2x_setiq(iq)は、出力値をスイッチしてくれるスイッチ関数である。初期化時点でデータをモジュールに溜め込んで、指示変数iqによりm_Igv2x_setiqでスイッチできる。メモリが足りなくなったらあきらめてファイルへ読み書きさせるメカにすることも可能であろう。モジュールを使う側からはそのメカを気にしなくて良いことになる。モジュール内サブルーチンの出力を引数で返すのは人間による整合性の二重記述が発生するのでよろしくない。データをモジュール内に貯め込むというコンセプトは重要であるー外から読めるが外からは書けない。
• スイッチ関数が不要なばあいも多い。on the flyでデータ生成してモジュール変数に戻すのもあり得る。

module m_igv2x !Return G vectors (integer sets) for given q points specifiec by qplist(:,iq)
  use m_struc_def,only: s_nv2
  public:: m_igv2xall_init, m_igv2x_setiq
  integer,protected,pointer,public :: igv2x(:,:)
  integer,protected,pointer,public :: napw,ndimh,ndimhx
  integer,allocatable,target,protected,public:: ndimhall(:)
  private
  integer,protected,private,target ::napw_z,ndimh_z,ndimhx_z
  type(s_nv2),allocatable,target,protected,private:: igv2xall(:)
  integer,allocatable,target,protected,private:: napwall(:),ndimhxall(:),igv2x_z(:,:)
  logical,private:: init=.true.
contains
  subroutine m_Igv2x_setiq(iq) ! Return G vectors for given qplist(:,iq)
    integer:: iq
    napw  => napwall(iq)
    ndimh => ndimhall(iq)  !nlmto+napw
    ndimhx=> ndimhxall(iq) !nlmto+napw (but x2 when SO=1)
    igv2x => igv2xall(iq)%nv2
  end subroutine m_Igv2x_setiq
  subroutine m_igv2xall_init(iqini,iqend) !initialization for qplist(iqini:iqend)  
    use m_qplist,only: qplist 
    use m_MPItk,only: master_mpi,procid,master
    integer:: iqini,iqend,iq
...

• またmoduleには引数なしでサブルーチンをcontainするやり方もある。これにはインターフェイス文の自動作成(.mod)という意味がある。変なネット資料だといちいち手でinterfaceを書けと指導しているものもある。複数のサブルーチン場合、何がpublicかを示すという意味もある。引数なしにすることにこだわらないといけない。変なmodule引数があるとどこで読み書きしているかわからなくなる。引数の整合性チェックをある程度してくれるのでとても助かる。optional変数とかname対応をさせたい場合は.modを作るのは必須。

• これ以外のmoduleの使い方としては「構造体格納module」というのもひとつ用いている。

逆行列計算のコード

moduleにいれてuseするのがいいのかもしれない。そのときにはmatinv(a)でよい。

subroutine matinv(n,a)
  implicit none
  integer :: n, info, ipiv(n)
  real(8):: a(n,n)
  real(8):: work(n*n)
  call dgetrf(n,n,a,n,ipiv,info)
  if(info/=0) then
     print *,' matinv: degtrf info=',info
     call exit(-1) !rx( ' matinv: degtrf ')
  endif
  call dgetri(n,a,n,ipiv,work,n*n,info)
  deallocate(work)
  if(info/=0) then
     print *,'matinv: degtri info=',info
     call exit(-1) !rx( 'matinv: degtri ')
  endif
end subroutine matinv

（rxはプリントしてexit(-1)する自作関数）

逆行列計算のコードその２

以下のような３点内挿を書くとき(wgt3が内挿重みを求める）、３ｘ３逆行列関数inverse33
が必要になった。

              do concurrent(it=itini:itend)     ! nt0p corresponds to efp
　　　　　　　　　　...
                 associate(x=>we_(it,itp),xi=>freq_r(ixs-1:ixs+1))!x=>we_ is \omega_\epsilon in Eq.(55). 
                   amat(1:3,1)= 1d0                    
                   amat(1:3,2)= xi(1:3)**2
                   amat(1:3,3)= xi(1:3)**4
                   wgt3(:,it,itp)= matmul([1d0,x**2,x**4], inverse33(amat)) 
                 endassociate

というわけで以下のようにしてみた。

  pure function inverse33(matrix) result(inverse) !Inverse of 3X3 matrix
    implicit none
    real(8),intent(in) :: matrix(3,3)
    real(8) :: inverse(3,3), det
    inverse(:,1)= crossf(matrix(:,2),matrix(:,3))
    inverse(:,2)= crossf(matrix(:,3),matrix(:,1))
    inverse(:,3)= crossf(matrix(:,1),matrix(:,2))
    det = sum(matrix(:,1)*inverse(:,1)) 
    inverse = transpose(inverse)
    inverse = inverse/det
  end function inverse33
  pure function crossf(a,b) result(c)
    implicit none
    intent(in):: a,b
    real(8):: a(3),b(3),c(3)
    c(1)=a(2)*b(3)-a(3)*b(2)
    c(2)=a(3)*b(1)-a(1)*b(3)
    c(3)=a(1)*b(2)-a(2)*b(1)
  end function crossf

展開しておいたほうがいいかもしれない。上述の逆行列コードもこういう形式のほうがいい。

mpiio

mpiioを用いて、以下のようなレコード長固定のrandom accessに似た形式にして利用している。

module m_mpiio !MPI-IO only for complex(8). Fixed length recl
  use m_nvfortran
  use mpi
  implicit none
  public:: openm,writem,readm,closem
  private
  integer,parameter::nfmax=100, nsize=16 !maxsize of opened file by openm
  integer:: ierr,fhl(nfmax)=-9999,iff=0
  integer(kind=mpi_offset_kind)::recll(nfmax)
contains
  function openm(newunit,file,recl) result(i) !recl=16*size
    integer::    newunit,     recl,info,amode
    character(*)::     file
    integer:: i
    info = mpi_info_null
    call mpi_file_open(MPI_COMM_WORLD, trim(file), mpi_mode_rdwr + mpi_mode_create,MPI_INFO_NULL, newunit,ierr)
    iff=iff+1
    if(iff>nfmax) call rx('m_mpiio:iff>nfmax')
    fhl(iff)   = newunit
    recll(iff) = recl !in byte
    i=0
  end function openm
  function writem(unit,rec,data) result(i)
    integer::unit
    integer(mpi_offset_kind) :: offset
    integer::rec,count
    complex(8):: data(1)
    integer:: i,ifx
    integer:: status(MPI_Status_size)
    ifx = fhl(findloc(unit==fhl,dim=1,value=.True.))
    offset= (rec-1)*recll(ifx)
    count = recll(ifx)/nsize     !    write(6,*)'writemmmmm',ifx,offset,count 
    call mpi_file_write_at(ifx, offset, data, count, MPI_DOUBLE_COMPLEX, status, ierr)
    i=0
  end function writem
  function readm(unit,rec,data) result(i)
    integer::unit
    integer::rec,count
    integer(mpi_offset_kind) :: offset
    integer::  status(MPI_STATUS_SIZE)
    complex(8):: data(1)
    integer:: i,ifx
    ifx = fhl(findloc(unit==fhl,dim=1,value=.True.))
    offset=(rec-1)*recll(ifx)
    count=recll(ifx)/nsize
    call mpi_file_read_at(ifx, offset, data, count, MPI_DOUBLE_COMPLEX, status, ierr)
    i=0
  end function readm
  function closem(unit) result(i)
    integer::unit
    integer:: i
    call mpi_file_close(unit, ierr)
    i=0
  end function closem
end module m_mpiio

compile option

たとえば、自分の場合、m_qplist.f90がコンパイラバグでコケるので、-O0でコンパイルするように指示している。コンパイラ依存をcppなどでソースに書き込むのは長期的に見るとホントしんどくなる。やめたほうがいい。

ifortの場合

やはりNaN初期化が安全側だと思う。

### for gfortran ###
mpif90 -Wsurprising -Wmaybe-uninitialized -O2  -g -fimplicit-none -finit-integer=NaN -finit-real=NaN -JOBJ.gfortran -IOBJ.gfortran -c ../main/qg4gw.m.f90 -o OBJ.gfortran/qg4gw.m.o

開発環境

• fortranでもコード開発の大半は手元のノートパソコンでやれる。フロントエンドを共用する大型システムは開発に向かないことも多い。ぼくの場合、ノートパソコン+ubuntu+gfortranでかなりの部分を行っている。最終仕上げを大型へもっていけばいい。emacs,vscodeを併用している。
• ワイドモニタ（A4で４枚分の幅がある）を2画面分割でやるのがよい。中央で割って200文字程度で2画面使える。-ffree-line-length-noneなどとして文字数制限は撤去している。コードの右側に！をつけてコメントを書くとみやすい。ノート＋HDMIのワイドモニタ。３モニタは使いきれんかった。

emacs

検索、置換、などに加えて以下の活用がおすすめ。知らない場合は試す.
ただ若い人ならemacsを試すより、vscodeが良いかもしれない。

• etag
• シェルウインドウのrename,別フレーム。複数のシェルを立ち上げる（ぼくの場合、必ずt1,t2とリネームしている）。
• gitモード　差分など
• ediffモード　差分と合成。
• grepモード
• コンパイルモードだとM-x compileと打ってコマンドを打つ。たとえば、

Compile command: cd ~/ecalj/SRC/exec;make -j lmf lmfa lmchk

などと打つ。タブで次の行。リターンでエラーへジャンプ。C-x b *compileでコンパイル画面に戻ってgで再コンパイル。shift+矢印で記憶が戻る.もちろん、c-x b c tabでcompileに戻れる。

• 「複数ファイルの一括変換」なども可能。
• カスタマイズに凝りすぎない。デフォルトを活用。

git,gitk

過去のソースc-x v l,リターン. dで差分表示。
差分のマージにはediffを使っている。
過去のソースcommitを指定してチェックアウト.
gitk --allも有効.