printf関数によって値が変化してしまう問題について

[aki]

関数f(x)=((Ax-y))^2のx[1], x[2]における勾配を求めたくて以下のようなコードを書きました。（ Aはn by nの行列、x,yはn次元ベクトル、 ((v))はベクトルvのノームを表しています。xでのiはxの添字です。）

コード:

#include <stdio.h>
#define dx 1e-7

double findGrad(double *a, double *y, double *x, int n, int i, double fm);
double func(double *a, double *y, double *x, int n);
void multiply(double *mat, double *x, int n);

int main(void){
	int n, i;
	n = 2;
	double a[n*n], x[n], y[n], df[n], fm;

	x[0] = 0; x[1] = 0;
	y[0] = 8; y[1] = 18;

	for(i=0;i<n*n;i++){
		a[i] = i+1;
	}

	fm = func(a,y,x,n);
	printf("%f\n", fm);

	for(i=0; i<n; i++){
			df[i] = findGrad(a, y, x, n, i, fm);
			printf("%f ", df[i]);
		}
	printf("\n");

}



double findGrad(double *a, double *y, double *x, int n, int i, double fm){
	double fx, xNew[n];
	int k;
	for(k=0; k<n; k++){
		xNew[k] = x[k];
		//printf("xNew: %f, %f\n", xNew[0], xNew[1]);
	}
	xNew[i] = x[i] + dx;
	//printf("xNew: %f\n", xNew[i]);
	fx = 0;
	//printf("fx: %f\n", fx);
	fx = func(a,y,xNew,n);
	//printf("fx: %f\n", (fx-fm)/dx);
	return (fx-fm)/dx;
}

double func(double *a, double *y, double *x, int n){
	int i;
	double v[n], norm = 0;

	multiply(a, x, n); //calculate a*x and put the result in x;
	for(i=0; i<n; i++){
		v[i] = x[i]-y[i];
		norm += v[i]*v[i];
	}
	return norm;
}

void multiply(double *mat, double *x, int n){
	double v[n], w[n];
	int i, j;

	for(i=0; i<n; i++){
		for(j=0; j<n; j++){
			v[j] = mat[i*n+j]*x[j];
			w[i] += v[j];
		}
	}
	//copy v to x
	for(i=0; i<n; i++){
		x[i] = w[i];
	}
}

チェック用に、n = 2, A = [1, 2; 3, 4], x = [0 ; 0], y = [8 ; 18]としたのでf(x)は以下の関数になります。
f(x) = 10x[1]^2 - 124x[1] +28x[1]x[2] + 20x[2]^2 - 176x[2] + 388
勾配は
df/dx[1] = 20x[1] + 28x[2] -124
df/dx[2] = 40x[2] + 28x[1] - 176　
であり、x[1] = 0, x[2] = 0のときそれぞれ-124, -176となるはずですが、実際には
x[1] = -123.999999
x[2] = -299.999994
となってしまい、x[2]のほうがなぜかx[1] + x[2]のように上書きされてしまっているようです。

しかしdouble findGrad(double *a, double *y, double *x, int n, int i, double fm)内各変数の値をprintfで確認すると以下のように正常な値を出力します。
x[1] = -123.999999
x[2] = -175.999999

色々な変数をprintfでチェックしてみたのですが、どのケースもprintfを使ったときは正常な値を出すのですが、printfを使わない場合は上書きされた値が出力されました。

なぜこのようなことが起こるのでしょうか。
どなたか御教示いただけると幸いです。

[超初級者]

ぱっと見、
multiply関数で
w[]
を初期化していないのに
+=
しているのは
まずいだろうな、と思います。

[へにっくす]

掲示されたコードでは、そのままコンパイルできないですね（少なくともVisual Studio 2013 Expressではコンパイルエラー）。

とりあえず以下のように修正させてもらいました。

コード:

#include <stdio.h>
#define dx (1e-7)

double findGrad(double *a, double *y, double *x, int n, int i, double fm);
double func(double *a, double *y, double *x, int n);
void multiply(double *mat, double *x, int n);

int main(void){
	int i;
	const int n = 2;
	double a[n*n], x[n], y[n], df[n], fm;

	x[0] = 0; x[1] = 0;
	y[0] = 8; y[1] = 18;

	for (i = 0; i<n*n; i++){
		a[i] = i + 1;
	}

	fm = func(a, y, x, n);
	printf("%f\n", fm);

	for (i = 0; i<n; i++){
		df[i] = findGrad(a, y, x, n, i, fm);
		printf("%f ", df[i]);
	}
	printf("\n");
}



double findGrad(double *a, double *y, double *x, int n, int i, double fm){
	double fx, *xNew;
	int k;
	xNew = new double[n];
	for (k = 0; k<n; k++){
		xNew[k] = x[k];
		//printf("xNew: %f, %f\n", xNew[i]);
	}
	xNew[i] = x[i] + dx;
	//printf("xNew: %f\n", xNew[i]);
	fx = 0;
	//printf("fx: %f\n", fx);
	fx = func(a, y, xNew, n);
	//printf("fx: %f\n", (fx-fm)/dx);
	delete[] xNew;
	return (fx - fm) / dx;
}

double func(double *a, double *y, double *x, int n){
	int i;
	double *v, norm = 0;

	v = new double[n];

	for (i = 0; i < n; i++) {
		v[i] = 0;
	}
	multiply(a, x, n); //calculate a*x and put the result in x;
	for (i = 0; i<n; i++){
		v[i] = x[i] - y[i];
		norm += v[i] * v[i];
	}

	delete[] v;
	return norm;
}

void multiply(double *mat, double *x, int n){
	double *v, *w;
	int i, j;

	v = new double[n];
	w = new double[n];

	for (i = 0; i < n; i++) {
		v[i] = 0;
		w[i] = 0;
	}
	for (i = 0; i<n; i++){
		for (j = 0; j<n; j++){
			v[j] = mat[i*n + j] * x[j];
			w[i] += v[j];
		}
	}
	//copy v to x
	for (i = 0; i<n; i++){
		x[i] = w[i];
	}

	delete[] v;
	delete[] w;
}

変えたところは、２行目の#defineをかっこつけたのと、10行目は定数なのでconstを付けました。
あとそれぞれの関数で変数nを使って配列を作成するところ(もちろん使った後はdeleteします)と、初期値を入れるようにしたことですね。
実行結果は以下の通り。printfあろうがなかろうが同じ結果になりましたよ。

printfなし

コード:

> ***.exe
388.000000
-123.999999 -175.999999
>

printfあり

コード:

> ***.exe
388.000000
xNew: 0.000000, 0.000000
xNew: 0.000000, 0.000000
xNew: 0.000000
fx: 0.000000
fx: -123.999999
-123.999999 xNew: 0.000000, 0.000000
xNew: 0.000000, 0.000000
xNew: 0.000000
fx: 0.000000
fx: -175.999999
-175.999999
>