C++で2次元の3次スプライン補間．方程式の解が不定の場合の対処法について

Question

C++で2次元の3次スプライン補間を行おうと考えています．
元々のcsvはこのBx.csv(　https://drive.google.com/file/d/1vcpLAAZIc6EL4QZ8YTfYthSwZu7UqGWR/view?usp=sharing　)というファイルを使用します．
このファイルの1行目がx座標，1列目がy座標を表しています．
それを2乗3乗した2次元配列を作成しそれを基に16本の連立方程式を作り，それを解くことによって0.01×0.01の各マス内における3次スプライン補間式の係数を導出したいです．
添付プログラム序盤のcsvの読み取りはこちら(　https://teratail.com/questions/289856　)を，連立方程式の解き方はこちら(　https://qiita.com/kubo_programmer/items/c42ec88bdba3c5c0919f　)を参考にさせていただきました．

これで計算したところcan not calc（自分で作成したエラー文です）と返されてしまいます．
おそらくcsvの中に値が十分に散らばっていない部分があり連立方程式の解が一位に定まらずエラーになってしまうのだと考えられます．
当然すべてのセルにとりあえず１を足すなどの処理は無意味で根本的な対策が必要みたいです．

このような場合に良い対処法などはありますでしょうか？

結果(2次元配列a～p)はこの後に別のプログラムで使いたいと思っていて，この計算がうまくいったらそれを別のcsvにでも格納して使いたいと思ってるので，最後の部分はこのままでいきたいです．

1日考えましたが現在の私の力ではどうにも解決できなさそうなので解決法をご教示いただけますと幸いです．
よろしくお願いします．

該当のソースコード

#include <iostream>
#include <string> 
#include <stdint.h>
#include <vector>
#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#include <iostream>
#include <math.h>
#include <fstream>
#include <sstream>
#include <algorithm> 
#include <chrono>
#include <utility> 
#include <time.h> 
#include <tuple> 
#include <cstdint> 
#include <cstdio>
//別のプログラムで使いたいのでincludeは極力このままで
using namespace std;


const int w = 400;
using Mat = vector<vector<double>>;
vector<double> parse_line(const string &line, char delim = ',')
{
	vector<double> values;

	istringstream ss(line);
	string token;
	while (getline(ss, token, delim))
		values.push_back(stod(token));
	return values;
}
Mat parse_csv(const string &path)
{
	ifstream ifs(path);
	string line;
	Mat mat;

	while (getline(ifs, line)) {
		// 1行分解析する。
		vector<double> values = parse_line(line);
		mat.push_back(values);
	}
	return mat;
}
double sq[w][w];
double cb[w][w];
double a[w][w], b[w][w], c[w][w], d[w][w], e[w][w], f[w][w], g[w][w], h[w][w];
double s[w][w], t[w][w], k[w][w], l[w][w], m[w][w], n[w][w], o[w][w], p[w][w];

int main() {
	

	Mat si;
	try {
		si = parse_csv("Bx.csv");
	}
	catch (exception &e) {
		cout << "Error: " << e.what() << "\n"; return 1;
	}

	int y=si.size();  // (縦の要素数)364
	int x= si.at(0).size(); //(横の要素数)184
	
	
	for(int i=0;i<y;i++) {
		for (int j = 0; j<x; j++) {
			sq[i][j] = si[i][j] * si[i][j];
			cb[i][j]= si[i][j] * si[i][j] * si[i][j];
		}
	}

	//ここまではうまくいっています．


	//f=a*xxx+b*xxxy+c*xxxyy+d*xxxyyy
	//+e*xx+f*xxy+g*xxyy+h*xxyyy
	//+s*x+t*xy+k*xyy+l*xyyy
	//+m+n*y+o*yy+p*yyy

	//j-1,j,j+1,j+2
	//i-1,1,2,3,4
	//i,5,6,7,8
	//i+1,9,10,11,12
	//i+2,13,14,15,16

	//多分ここからが問題と思われます
	
	//1つずつ解く
	for (int i = 2; i <= y - 3; i++) {
		for (int j = 2; j <= x - 3; j++) {
			//このループを各２２～１００などとしても同様のエラーが出るので
			//問題はこのループの数字ではないと思います．
			
			//連立方程式を解く

			double z[16][17] = {
				{cb[0][j - 1],cb[0][j - 1] * si[i - 1][0],cb[0][j - 1] * sq[i - 1][0],cb[0][j - 1] * cb[i - 1][0],sq[0][j - 1], sq[0][j - 1] * si[i - 1][0], sq[0][j - 1] * sq[i - 1][0], sq[0][j - 1] * cb[i - 1][0], si[0][j - 1], si[0][j - 1] * si[i - 1][0], si[0][j - 1] * sq[i - 1][0], si[0][j - 1] * cb[i - 1][0], 1,  si[i - 1][0], sq[i - 1][0],  cb[i - 1][0],si[i-1][j-1]},//1
				{ cb[0][j],cb[0][j ] * si[i - 1][0],cb[0][j ] * sq[i - 1][0],cb[0][j ] * cb[i - 1][0],sq[0][j], sq[0][j ] * si[i - 1][0], sq[0][j ] * sq[i - 1][0], sq[0][j ] * cb[i - 1][0],si[0][j], si[0][j] * si[i - 1][0], si[0][j ] * sq[i - 1][0], si[0][j ] * cb[i - 1][0],1,  si[i - 1][0], sq[i - 1][0],  cb[i - 1][0], si[i - 1][j ] },//2
				{ cb[0][j + 1],cb[0][j + 1] * si[i - 1][0],cb[0][j + 1] * sq[i - 1][0],cb[0][j + 1] * cb[i - 1][0],sq[0][j + 1], sq[0][j + 1] * si[i - 1][0], sq[0][j + 1] * sq[i - 1][0], sq[0][j + 1] * cb[i - 1][0],si[0][j + 1], si[0][j + 1] * si[i - 1][0], si[0][j + 1] * sq[i - 1][0], si[0][j + 1] * cb[i - 1][0],1,si[i - 1][0], sq[i - 1][0],  cb[i - 1][0] , si[i - 1][j + 1] },//3
				{ cb[0][j + 2],cb[0][j + 2] * si[i - 1][0],cb[0][j + 2] * sq[i - 1][0],cb[0][j + 2] * cb[i - 1][0],sq[0][j + 2], sq[0][j + 2] * si[i - 1][0], sq[0][j + 2] * sq[i - 1][0], sq[0][j + 2] * cb[i - 1][0],si[0][j + 2], si[0][j + 2] * si[i - 1][0], si[0][j + 2] * sq[i - 1][0], si[0][j + 2] * cb[i - 1][0],1,  si[i - 1][0], sq[i - 1][0],  cb[i - 1][0], si[i - 1][j +2] },//4
				{ cb[0][j - 1],cb[0][j - 1] * si[i ][0],cb[0][j - 1] * sq[i ][0],cb[0][j - 1] * cb[i ][0],sq[0][j - 1], sq[0][j - 1] * si[i ][0], sq[0][j - 1] * sq[i ][0], sq[0][j - 1] * cb[i ][0],si[0][j - 1], si[0][j - 1] * si[i ][0], si[0][j - 1] * sq[i ][0], si[0][j - 1] * cb[i ][0],1,  si[i ][0], sq[i ][0],  cb[i ][0],si[i ][j - 1] },//5
				{ cb[0][j ],cb[0][j ] * si[i][0],cb[0][j ] * sq[i][0],cb[0][j ] * cb[i][0],	sq[0][j ], sq[0][j ] * si[i][0], sq[0][j] * sq[i][0], sq[0][j ] * cb[i][0],si[0][j ], si[0][j ] * si[i][0], si[0][j] * sq[i][0], si[0][j ] * cb[i][0],1,si[i][0], sq[i][0],  cb[i][0],si[i][j] },//6
				{ cb[0][j + 1],cb[0][j + 1] * si[i][0],cb[0][j + 1] * sq[i][0],cb[0][j + 1] * cb[i][0],sq[0][j +1], sq[0][j + 1] * si[i][0], sq[0][j +1] * sq[i][0], sq[0][j + 1] * cb[i][0],si[0][j + 1], si[0][j + 1] * si[i][0], si[0][j + 1] * sq[i][0], si[0][j + 1] * cb[i][0],1,si[i][0], sq[i][0],  cb[i][0],si[i][j + 1] },//7
				{ cb[0][j + 2],cb[0][j + 2] * si[i][0],cb[0][j + 2] * sq[i][0],cb[0][j + 2] * cb[i][0],sq[0][j + 2], sq[0][j + 2] * si[i][0], sq[0][j + 2] * sq[i][0], sq[0][j + 2] * cb[i][0],si[0][j + 2], si[0][j + 2] * si[i][0], si[0][j + 2] * sq[i][0], si[0][j + 2] * cb[i][0], 1,  si[i][0], sq[i][0],  cb[i][0],si[i][j + 2] },//8
				{ cb[0][j - 1],cb[0][j - 1] * si[i + 1][0],cb[0][j - 1] * sq[i + 1][0],cb[0][j - 1] * cb[i + 1][0],sq[0][j - 1], sq[0][j - 1] * si[i + 1][0], sq[0][j - 1] * sq[i + 1][0], sq[0][j - 1] * cb[i + 1][0],si[0][j - 1], si[0][j - 1] * si[i + 1][0], si[0][j - 1] * sq[i + 1][0], si[0][j - 1] * cb[i + 1][0],1,  si[i - 1][0], sq[i + 1][0],  cb[i + 1][0],si[i + 1][j - 1] },//9
				{ cb[0][j],cb[0][j] * si[i + 1][0],cb[0][j] * sq[i + 1][0],cb[0][j] * cb[i + 1][0],sq[0][j], sq[0][j] * si[i + 1][0], sq[0][j] * sq[i + 1][0], sq[0][j] * cb[i + 1][0],si[0][j], si[0][j] * si[i + 1][0], si[0][j] * sq[i + 1][0], si[0][j] * cb[i + 1][0],1,  si[i + 1][0], sq[i + 1][0],  cb[i + 1][0], si[i +1][j] },//10
				{ cb[0][j + 1],cb[0][j + 1] * si[i + 1][0],cb[0][j + 1] * sq[i + 1][0],cb[0][j + 1] * cb[i + 1][0],sq[0][j + 1], sq[0][j + 1] * si[i + 1][0], sq[0][j + 1] * sq[i + 1][0], sq[0][j + 1] * cb[i + 1][0],si[0][j + 1], si[0][j + 1] * si[i + 1][0], si[0][j + 1] * sq[i + 1][0], si[0][j + 1] * cb[i + 1][0],1,si[i + 1][0], sq[i + 1][0],  cb[i + 1][0] , si[i + 1][j + 1] },//11
				{ cb[0][j + 2],cb[0][j + 2] * si[i + 1][0],cb[0][j + 2] * sq[i + 1][0],cb[0][j + 2] * cb[i + 1][0],sq[0][j + 2], sq[0][j + 2] * si[i + 1][0], sq[0][j + 2] * sq[i + 1][0], sq[0][j + 2] * cb[i + 1][0],si[0][j + 2], si[0][j + 2] * si[i + 1][0], si[0][j + 2] * sq[i + 1][0], si[0][j + 2] * cb[i + 1][0],1,si[i + 1][0], sq[i+1][0],  cb[i + 1][0],si[i + 1][j + 2] },//12
				{ cb[0][j - 1],cb[0][j - 1] * si[i + 2][0],cb[0][j - 1] * sq[i + 2][0],cb[0][j - 1] * cb[i + 2][0],sq[0][j - 1], sq[0][j - 1] * si[i + 2][0], sq[0][j - 1] * sq[i + 2][0], sq[0][j - 1] * cb[i + 2][0],si[0][j - 1], si[0][j - 1] * si[i + 2][0], si[0][j - 1] * sq[i + 2][0], si[0][j - 1] * cb[i + 2][0],1, si[i +2][0], sq[i + 2][0],  cb[i + 2][0],si[i + 2][j - 1] },//13
				{ cb[0][j],cb[0][j] * si[i + 2][0],cb[0][j] * sq[i + 2][0],cb[0][j] * cb[i + 2][0],sq[0][j], sq[0][j] * si[i + 2][0], sq[0][j] * sq[i + 2][0], sq[0][j] * cb[i + 2][0], si[0][j], si[0][j] * si[i + 2][0], si[0][j] * sq[i + 2][0], si[0][j] * cb[i + 2][0], 1,  si[i + 2][0], sq[i + 2][0],  cb[i + 2][0], si[i + 2][j] },//14
				{ cb[0][j + 1],cb[0][j + 1] * si[i + 2][0],cb[0][j + 1] * sq[i + 2][0],cb[0][j + 1] * cb[i + 2][0],sq[0][j + 1], sq[0][j + 1] * si[i + 2][0], sq[0][j + 1] * sq[i + 2][0], sq[0][j + 1] * cb[i + 2][0],si[0][j + 1], si[0][j + 1] * si[i + 2][0], si[0][j + 1] * sq[i + 2][0], si[0][j + 1] * cb[i + 2][0],1, si[i + 2][0], sq[i + 2][0],  cb[i + 2][0] , si[i + 2][j + 1] },//15
				{ cb[0][j + 2],cb[0][j + 2] * si[i + 2][0],cb[0][j + 2] * sq[i + 2][0],cb[0][j + 2] * cb[i + 2][0],sq[0][j + 2], sq[0][j + 2] * si[i + 2][0], sq[0][j + 2] * sq[i + 2][0], sq[0][j + 2] * cb[i + 2][0],si[0][j + 2], si[0][j + 2] * si[i + 2][0], si[0][j + 2] * sq[i + 2][0], si[0][j + 2] * cb[i + 2][0], 1,  si[i + 2][0], sq[i + 2][0],  cb[i + 2][0],si[i + 2][j + 2] }//16
			};
			
			for (int ii = 0; ii<16; ii++) {
				int maxLine = ii;
				for (int jj = ii; jj<16; jj++) if (abs(z[maxLine][ii])<abs(z[jj][ii])) maxLine = jj;
				for (int jj = ii; jj<= 16; jj++) swap(z[ii][jj], z[maxLine][jj]);

				 double beginNum = z[ii][ii];
				if (!beginNum) {
					cout << "can not calc" << endl;
					return 0;
				}

				for (int jj = ii; jj <= 16; jj++)z[ii][jj] /= beginNum;

				for (int jj = ii + 1; jj<16; jj++) {
					 double beginDelete = z[jj][ii];
					for (int k = ii; k <= 16; k++)z[jj][k] -= beginDelete*z[ii][k];
				}
			}
			for (int ii = 15; ii >= 0; ii--) {
				for (int jj = ii + 1; jj<16; jj++)z[ii][16] -= z[jj][16] * z[ii][jj];
			}
			a[i][j] = z[0][16];
			b[i][j] = z[1][16];
			c[i][j] = z[2][16];
			d[i][j] = z[3][16];
			e[i][j] = z[4][16];
			f[i][j] = z[5][16];
			g[i][j] = z[6][16];
			h[i][j] = z[7][16];
			s[i][j] = z[8][16];
			t[i][j] = z[9][16];
			k[i][j] = z[10][16];
			l[i][j] = z[11][16];
			m[i][j] = z[12][16];
			n[i][j] = z[13][16];
			o[i][j] = z[14][16];
			p[i][j] = z[15][16];
			
		}
	}
	

	cout << o[244][55];//無事に解析できたかとりあえず確認



	//この問題が解決後自分でここに続きを打ち込む予定です

	return 0;
}

Comments

[cateye]

double a[w][w], b[w][w], c[w][w], d[w][w], e[w][w], f[w][w], g[w][w], h[w][w];
double s[w][w], t[w][w], k[w][w], l[w][w], m[w][w], n[w][w], o[w][w], p[w][w];
・・・で、何MBになりますか？・・・スタックにこれだけ取れるのでしょうか？

[Kinsho]

ご回答ありがとうございます．家庭用ゲーミングPCで弾かれるのは当然予想していましたが数値計算を専門にしている研究室の計算機でもエラーを返されてしまいました．
これくらいの容量ならスタックは問題なく置けるスペックだと思われます．

[mjk]

（研究室の計算機のスタックが特殊なのかは分かりませんが）
cateyeさんがご指摘されていますが、
double sq[w][w];
double cb[w][w];
その他配列をグローバル宣言にしてみるなどを試してみましたか？

私も似たような経験がありもしかしたらそれで解決するかもしれません。
https://teratail.com/questions/285169
配列サイズが大きくなりすぎるとメモリ不足になるのでそちらが原因では無いかとのご指摘かと思います。
それでも同様のエラーが出るようなら範囲外の配列を参照しているなどが疑われます。
配列[10]しかないのに配列[11]を参照しようとしたなど。

[mjk]

質問コードで試したところ結果は何も表示されませんでしたが配列をグローバル宣言したらセグメンテーションエラーは無くなりました。
https://gyazo.com/2205dd10dba524d3d4347fc522498df9

[Kinsho]

ご指摘ありがとうございます．ご指摘の通り配列宣言を全てグローバル変数の外に出したらその問題は解決できました．素人ゆえにPCのスペックを盲信していた私が間違っていました．

[mjk]

いえいえ初心者の私も直感で動くはずなのに何故動かないのか分からずに戸惑いました。
質問して回答してくださった方に教えて頂き初めて理解しました。
回答してくださった方に感謝しております。

Answer 1

皆様のご指摘の個所を修正後連立方程式の不定解問題の解消(厳密には近似)方法について自己解決ができましたのでのちに同じ問題に直面する方のために共有しておきます．

csvのxy座標に同じものが入っているとrank16を確保できないので，x=0,0.01,0.02,...のように始まる場合は最初の列をコピペした後にx=「-0.00000001」,0,0.01,0.02,...のようにするとrank16を確保できて16本の連立方程式で解が求まる．またこのデータの縦列のように中間に0をまたぐ場合はその0を同様に0.00000001などとしてやるといいと思います．

Answer 2

iやjを使ったforループの中に，iやjを使ったforが記述されているので，
その動作は想定通りではないものと推測されます．

Comments

[Kinsho]

確かにそこはうっかりしていました．ご指摘ありがとうございます．
ただi,jループ内にあったi,jをとりあえずii,jjなどとしてみましたが同様のエラーを返されてしまいました．．．