洞川温泉旅行記

今年の夏休みは奥さんと一緒に奈良県洞川温泉に行き、1泊しました。 去年の12月にはじめて行ったとき、空気が澄んでいてとても良い雰囲気でしたので、夏はどうかなと思い行ってみました。

まず初日は新幹線で京都まで行き、京都駅でレンタカーを借りました。 京都で櫛のお店、それから、お香と雑貨のお店を見ました。雨が少し降りはじめていましたが、歩いているとじっとりと汗をかく天気でした。 奥さんが赤い櫛とお線香を買えたのでよかったです。

その後、お腹が減ったので、お香と雑貨のお店の近くにある料理屋さんで軽く食べました。 僕はむぎとろご飯を食べ、奥さんは抹茶くず餅を食べました。 抹茶くず餅を少し分けてもらったのですが、もちゃもちゃとした食感でおいしかったです。

涼んでお腹がいっぱいになったところで、早速、奈良県洞川温泉の宿へ出発しました。 国道309号線の運転難易度が高いというようなことをネット記事で読んでいたので焦りました。 僕は運転が下手でよく慌ててしまいます。また、地図がよく読めません。 ナビに頼りきりなのですが、ナビの指示どおりで果たして到着できるのかが心配でした。 なんとかナビの支持どおりに運転して2時間半くらいで到着できました。 途中でトンネルを何個も通り、急カーブをなんとかこなしました。 真っ暗なトンネルを運転しているとき、ジージーッというようなたくさんの蝉の声のようなうなりを聞きました。 このあたりには心霊スポットのトンネルがあることを聞いていましたので少し怖くなりました。

そうこうしてなんとか宿に着き、チェックインしました。そのとき、おいしい赤いジュースをいただきました。 宿の夕食はヤマメの刺し身、天ぷら、そうめん、ハモなどおいしくいただきました。 奥さんの分もかなり食べたのでお腹が限界でした。

次の日は、宿の近くの川(山上川というらしい)でニジマスにエサをあげました。 川に入ったのですが、奥さんがずんずんと先へ行ってしまい、足の裏が痛くて歩くのが大変でした。 川底の石には苔が生えており、すべりそうでした。ニジマスはこの苔を食べているのでしょうか。 でも、坂の下のお店で売っていたニジマスのエサを橋の上から撒くと、鯉のように大量のニジマスが出てきたので、 このエサだけで十分かもしれません。

その後、去年に続き、今年も龍泉寺にお参りし、京都に戻りました。 祇園の豊田愛山堂というお店で、竜涎香の匂袋を買いました。 竜涎香の甘いような不思議な匂いが好きになりました。

引き戻しと押し出し

多様体間の写像によって誘導される引き戻しと押し出しという写像に対して、以下の関係が成り立つことを示したい。

1. 関数の引き戻し

多様体 Mから多様体 Nへの写像 fを考える。

 g\colon N \to \mathbb{R}に対して

が定義できる。これを g fによる M上への引き戻しという。

2. ベクトル場の積分曲線

多様体  M (m次元)上のベクトル場

に沿った曲線を考える。

そのような曲線をベクトル場の積分曲線といい、以下で定義される。

曲線 c\colon \mathbb{R} \to M

がベクトル場

積分曲線または解曲線であるとは、

が成り立つことである。

ベクトル場の積分曲線は重要である。なぜなら、

力学の問題は, どの定式化にせよ, 基礎方程式によって与えられたベクトル場にたいする積分曲線を求めるという形で表せる 山本義隆, 中村孔一 (1998) 「解析力学1 (朝倉物理学大系)」

3. ベクトル場の押し出し

 p = c(\tau) \in M とする。  f\colon M \to N微分同相写像のとき、写像

 f によるベクトル場  X の押し出し、または、  f微分写像という。 ベクトル場は微分作用素でもあるので  g\colon N \to \mathbb{R} への作用を考える。

よって、以下が成り立つ。

ラプラス・ベルトラミ作用素を使った3次元極座標のラプラシアンの導出

前回の記事では曲がった空間のラプラシアンであるラプラス・ベルトラミ作用素を求めた。

nyonyo.hatenablog.com

今回は以下のラプラス・ベルトラミ作用素を使って3次元極座標ラプラシアンを導出してみる。

これを3次元極座標で計算するには3次元極座標の計量を求めなければならない。 以下のように線素を使って求めることができる。

この計量をラプラス・ベルトラミ作用素に代入して計算すれば、以下のように3次元極座標ラプラシアンが導出できる。

ラプラス・ベルトラミ作用素(ラプラシアンの一般化)

曲がった空間のラプラシアンであるラプラス・ベルトラミ作用素を求める。 これを使うと3次元極座標でのラプラシアンが簡単に計算できる。

ラプラス・ベルトラミ作用素は以下。

ポイントはラプラシアンに含まれる微分を以下の共変微分で置き換えること。

置き換えて、以下のようになる。

(5)→(6)の式変形について説明する。 まず、共変微分に含まれるクリストッフェル記号は次のように計量を使って書ける。 ラプラシアンでは上付き添字と下付き添字が同じである(10)〜(13)を使う。

(12)→(13)のように式変形できる理由を以下に示す。

次回はこのラプラス・ベルトラミ作用素を使って実際に3次元極座標でのラプラシアンの表示を導きたい。

ARMv8.3-CompNum FCMLA命令の動作

FCMLA命令の動作について書く。

FCMLA <Zda>.<T>, <Pg>/M, <Zn>.<T>, <Zm>.<T>, <const>

複素数の乗算はFCMLA命令2つで書ける。 ベクトルレジスタには複数の複素数を格納できる。 偶数番目が実部で、奇数番目が虚部。 ベクトルレジスタが512ビットの場合、倍精度(double, 64ビット)の複素数を4つ格納できる。

  • 複素数1: 実部が0番目, 虚部が1番目
  • 複素数2: 実部が2番目, 虚部が3番目
  • 複素数3: 実部が4番目, 虚部が5番目
  • 複素数4: 実部が6番目, 虚部が7番目

この前提のもと、FCMLA命令の動作は以下の表のようになる。

const 実部 虚部
0 Re(Zda)+=Re(Zn)*Re(Zm) Im(Zda)+=Re(Zn)*Im(Zm)
90 Re(Zda)+=Im(Zn)*(-Im(Zm)) Im(Zda)+=Im(Zn)*Re(Zm)
180 Re(Zda)+=Re(Zn)*(-Re(Zm)) Im(Zda)+=Re(Zn)*(-Im(Zm))
270 Re(Zda)+=Im(Zn)*Im(Zm) Im(Zda)+=Im(Zn)*(-Re(Zm))

確認してみる。

#include <stdio.h>
#include <arm_sve.h>

int main() {
  double a[8];
  double b[8];
  double c[8];

  for (int i=0; i<8; i++)
    a[i] = i % 2 == 0 ? -i : i;

  for (int i=0; i<8; i++)
    b[i] = i * 2.0;

  for (int i=0; i<8; i++)
    c[i] = 0.0;

  svbool_t pg = svptrue_b64();

  svfloat64_t avec = svld1(pg, a);
  svfloat64_t bvec = svld1(pg, b);
  svfloat64_t cvec = svld1(pg, c);

  puts("*** rotate == 0 ***");
  svfloat64_t result = svcmla_x(pg, cvec, avec, bvec, 0);
  svst1(pg, c, result);
  for (int i=0; i<8; i++)
    printf("c[%d]=%f\n", i, c[i]);

  puts("check:");
  for (int i=0; i<8; i+=2) {
    // re
    c[i] = a[i] * b[i];
    // im
    c[i+1] = a[i] * b[i+1];
  }
  for (int i=0; i<8; i++)
    printf("c[%d]=%f\n", i, c[i]);

  puts("*** rotate == 90 ***");
  result = svcmla_x(pg, cvec, avec, bvec, 90);
  svst1(pg, c, result);
  for (int i=0; i<8; i++)
    printf("c[%d]=%f\n", i, c[i]);

  puts("check:");
  for (int i=0; i<8; i+=2) {
    // re
    c[i] = a[i+1] * -b[i+1];
    // im
    c[i+1] = a[i+1] * b[i];
  }
  for (int i=0; i<8; i++)
    printf("c[%d]=%f\n", i, c[i]);

  puts("*** rotate == 180 ***");
  result = svcmla_x(pg, cvec, avec, bvec, 180);
  svst1(pg, c, result);
  for (int i=0; i<8; i++)
    printf("c[%d]=%f\n", i, c[i]);

  puts("check:");
  for (int i=0; i<8; i+=2) {
    // re
    c[i] = a[i] * -b[i];
    // im
    c[i+1] = a[i] * -b[i+1];
  }
  for (int i=0; i<8; i++)
    printf("c[%d]=%f\n", i, c[i]);

  puts("*** rotate == 270 ***");
  result = svcmla_x(pg, cvec, avec, bvec, 270);
  svst1(pg, c, result);
  for (int i=0; i<8; i++)
    printf("c[%d]=%f\n", i, c[i]);

  puts("check:");
  for (int i=0; i<8; i+=2) {
    // re
    c[i] = a[i+1] * b[i+1];
    // im
    c[i+1] = a[i+1] * -b[i];
  }
  for (int i=0; i<8; i++)
    printf("c[%d]=%f\n", i, c[i]);
}
$ aarch64-linux-gnu-gcc-10 fcmla.c -march=armv8.3-a+sve -static
$ qemu-aarch64 -cpu max,sve-max-vq=4 a.out
*** rotate == 0 ***
c[0]=0.000000
c[1]=0.000000
c[2]=-8.000000
c[3]=-12.000000
c[4]=-32.000000
c[5]=-40.000000
c[6]=-72.000000
c[7]=-84.000000
check:
c[0]=0.000000
c[1]=0.000000
c[2]=-8.000000
c[3]=-12.000000
c[4]=-32.000000
c[5]=-40.000000
c[6]=-72.000000
c[7]=-84.000000
*** rotate == 90 ***
c[0]=-2.000000
c[1]=0.000000
c[2]=-18.000000
c[3]=12.000000
c[4]=-50.000000
c[5]=40.000000
c[6]=-98.000000
c[7]=84.000000
check:
c[0]=-2.000000
c[1]=0.000000
c[2]=-18.000000
c[3]=12.000000
c[4]=-50.000000
c[5]=40.000000
c[6]=-98.000000
c[7]=84.000000
*** rotate == 180 ***
c[0]=0.000000
c[1]=0.000000
c[2]=8.000000
c[3]=12.000000
c[4]=32.000000
c[5]=40.000000
c[6]=72.000000
c[7]=84.000000
check:
c[0]=-0.000000
c[1]=-0.000000
c[2]=8.000000
c[3]=12.000000
c[4]=32.000000
c[5]=40.000000
c[6]=72.000000
c[7]=84.000000
*** rotate == 270 ***
c[0]=2.000000
c[1]=0.000000
c[2]=18.000000
c[3]=-12.000000
c[4]=50.000000
c[5]=-40.000000
c[6]=98.000000
c[7]=-84.000000
check:
c[0]=2.000000
c[1]=-0.000000
c[2]=18.000000
c[3]=-12.000000
c[4]=50.000000
c[5]=-40.000000
c[6]=98.000000
c[7]=-84.000000

0.0-0.0でずれているところがあるのはなぜ?

C++でrestrictと同様の機能を実現する方法

ループのベクトル化や、その他の最適化を促進させることを目的として、C言語のrestrict修飾子を使うことがよくある。

restrict修飾子を指定したポインタは、同じ関数の他のポインタが指すメモリアドレスを指さない。 複数のポインタが同じメモリアドレスを指すと、メモリ依存関係によって、ベクトル化ができなくなることが多い。 restrictを使うことで、コンパイラが「メモリ依存が無い」と判定できるようになって、ベクトル化が可能になることがある。

以下のようにするとベクトル化されない。 コンパイラはruntimeチェックを入れることで、restrictが無いプログラムもベクトル化してしまうので、-mllvm -pragma-vectorize-memory-check-threshold=0を指定している。

norestrict.c

void test(double *a, double *b, int len) {
  for (int i=1; i<len; i++) {
    a[i] *= b[i] + i;
  }
}
$ clang -O2 -S norestrict.c -fno-unroll-loops -Rpass=loop-vectorize -mavx512f -mllvm -pragma-vectorize-memory-check-threshold=0

以下のようにrestrictを付けるとベクトル化されることがメッセージからわかる。

restrict.c

void test(double *restrict a, double *restrict b, int len) {
  for (int i=1; i<len; i++) {
    a[i] *= b[i] + i;
  }
}
$ clang -O2 -S restrict.c -fno-unroll-loops -Rpass=loop-vectorize -mavx512f -mllvm -pragma-vectorize-memory-check-threshold=0
restrict.c:2:3: remark: vectorized loop (vectorization width: 8, interleaved count: 1) [-Rpass=loop-vectorize]
  for (int i=1; i<len; i++) {
  ^

一方、C++の場合は、以下のように構造体のメンバとすることでベクトル化できるようになる。 これは、「別のタグ名を持つ構造体への2つのポインタは別名ではない」というルールを使っている。

a.cpp

struct A {
  double v;
};

struct B {
  double v;
};

void test(struct A *as, struct B *bs, int len) {
  for (int i=1; i<len; i++) {
    as[i].v *= bs[i].v + i;
  }
}
$ clang++ -O2 -S a.cpp -fno-unroll-loops -Rpass=loop-vectorize -mavx512f -mllvm -pragma-vectorize-memory-check-threshold=0
a.cpp:10:3: remark: vectorized loop (vectorization width: 8, interleaved count: 1) [-Rpass=loop-vectorize]
  for (int i=1; i<len; i++) {
  ^

上は、-fno-strict-aliasingを指定すると、ベクトル化されなくなる。

$ clang++ -O2 -S a.cpp -fno-unroll-loops -Rpass=loop-vectorize -mavx512f -mllvm -pragma-vectorize-memory-check-threshold=0 -fno-strict-aliasing

また、C++では、コンパイラの独自拡張である__restrict____restrictが使える様子。

restrict.cpp

void test(double *__restrict__ a, double *__restrict__ b, int len) {
  for (int i=1; i<len; i++) {
    a[i] *= b[i] + i;
  }
}
 clang++ -O2 -S restrict.cpp -fno-unroll-loops -Rpass=loop-vectorize -mavx512f -mllvm -pragma-vectorize-memory-check-threshold=0 
restrict.cpp:2:3: remark: vectorized loop (vectorization width: 8, interleaved count: 1) [-Rpass=loop-vectorize]
  for (int i=1; i<len; i++) {
  ^

低レベルプログラミング | Igor Zhirkov, 吉川 邦夫, 吉川 邦夫 |本 | 通販 | Amazonのp.351 「14.6 厳密な別名のルール」を参考にした。