[cig-commits] [commit] QA: Use proper C indexing. (37f4095)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon Jan 20 11:49:31 PST 2014
Repository : ssh://geoshell/specfem2d
On branch : QA
Link : https://github.com/geodynamics/specfem2d/compare/28743f19b9f9fdb75d359c135053825c0ffd05b3...5e8aa55e68fd17b6f475fb65531b84195e497aa1
>---------------------------------------------------------------
commit 37f40954e5da6b5069256faea98547a02edf5dab
Author: Elliott Sales de Andrade <esalesde at physics.utoronto.ca>
Date: Mon Jan 6 03:06:57 2014 -0500
Use proper C indexing.
No more of this pseudo-Fortran indexing. It's odd and causes unnecessary
work (to be fixed).
>---------------------------------------------------------------
37f40954e5da6b5069256faea98547a02edf5dab
src/specfem2D/attenuation_compute_param.c | 250 +++++++++++++++---------------
1 file changed, 124 insertions(+), 126 deletions(-)
diff --git a/src/specfem2D/attenuation_compute_param.c b/src/specfem2D/attenuation_compute_param.c
index a16784d..d77a1aa 100644
--- a/src/specfem2D/attenuation_compute_param.c
+++ b/src/specfem2D/attenuation_compute_param.c
@@ -16,10 +16,10 @@
void constant_Q2_sub(double f1, double f2, int n, double Q, double *tau_s, double *tau_e);
void nrerror(const char *error_text);
-double *dvector(int nl, int nh);
-double **dmatrix(int nrl, int nrh, int ncl, int nch);
-void free_dvector(double *v, int nl, int nh);
-void free_dmatrix(double **m, int nrl, int nrh, int ncl, int nch);
+double *dvector(int n);
+double **dmatrix(int nr, int nc);
+void free_dvector(double *v, int n);
+void free_dmatrix(double **m, int nr, int nc);
void dsvdcmp(double **a, int m, int n, double *w, double **v);
void initialize(double f1, double f2, int n, double Q, double *x1, double *x2);
void derivatives(double f1, double f2, int n, double Q, double *x1, double *x2, double *gradient, double **hessian);
@@ -66,46 +66,46 @@ FC_FUNC_(attenuation_compute_param,ATTENUATION_COMPUTE_PARAM)(int *nmech_in,
exit(1);
}
-/* loop on the Q1 dilatation mode (nu = 1) and Q2 shear mode (nu = 2) defined in Carcione's papers */
- for (nu = 1; nu <= 2; nu++) {
+/* loop on the Q1 dilatation mode (nu = 0) and Q2 shear mode (nu = 1) defined in Carcione's papers */
+ for (nu = 0; nu < 2; nu++) {
/* assign Q1 or Q2 to generic variable Q_value which is used for the calculations */
- if (nu == 1) { Q_value = target_Q1 ; }
- if (nu == 2) { Q_value = target_Q2 ; }
+ if (nu == 0) { Q_value = target_Q1 ; }
+ if (nu == 1) { Q_value = target_Q2 ; }
/* no need to compute these parameters if there is no attenuation; it could lead to a division by zero in the code */
if (Q_value > 0.00001) {
- tau_s = dvector(1, n);
- tau_e = dvector(1, n);
+ tau_s = dvector(n);
+ tau_e = dvector(n);
constant_Q2_sub(f1, f2, n, Q_value, tau_s, tau_e);
/* output in Fortran90 format */
- for (i = 1; i <= n; i++) {
+ for (i = 0; i < n; i++) {
/* We put the results in tau_sigma_nu to get them in fortran. */
- if ( nu == 1 ) {
- tau_sigma_nu1[i-1] = tau_s[i];
+ if ( nu == 0 ) {
+ tau_sigma_nu1[i] = tau_s[i];
}
- if ( nu == 2 ) {
- tau_sigma_nu2[i-1] = tau_s[i];
+ if ( nu == 1 ) {
+ tau_sigma_nu2[i] = tau_s[i];
}
}
- for (i = 1; i <= n; i++) {
+ for (i = 0; i < n; i++) {
/* We put the results in tau_epsilon_nu to get them in fortran. */
- if ( nu == 1 ) {
- tau_epsilon_nu1[i-1] = tau_e[i];
+ if ( nu == 0 ) {
+ tau_epsilon_nu1[i] = tau_e[i];
}
- if ( nu == 2 ) {
- tau_epsilon_nu2[i-1] = tau_e[i];
+ if ( nu == 1 ) {
+ tau_epsilon_nu2[i] = tau_e[i];
}
}
- free_dvector(tau_s, 1, n);
- free_dvector(tau_e, 1, n);
+ free_dvector(tau_s, n);
+ free_dvector(tau_e, n);
}
}
}
@@ -129,15 +129,15 @@ void constant_Q2_sub(double f1, double f2, int n, double Q, double *tau_s, doubl
exit(1);
}
- x1 = dvector(1, n);
- x2 = dvector(1, n);
- gradient = dvector(1, n);
- hessian = dmatrix(1, n, 1, n);
- for(i=1;i<=n;i++) {
+ x1 = dvector(n);
+ x2 = dvector(n);
+ gradient = dvector(n);
+ hessian = dmatrix(n, n);
+ for(i=0;i<n;i++) {
x1[i]=0.0;
x2[i]=0.0;
gradient[i]=0.0;
- for(j=1;j<=n;j++) hessian[i][j]=0.0;
+ for(j=0;j<n;j++) hessian[i][j]=0.0;
}
initialize(f1, f2, n, Q, x1, x2);
@@ -146,18 +146,18 @@ void constant_Q2_sub(double f1, double f2, int n, double Q, double *tau_s, doubl
invert(x1, gradient, hessian, n);
- free_dvector(gradient, 1, n);
- free_dmatrix(hessian, 1, n, 1, n);
+ free_dvector(gradient, n);
+ free_dmatrix(hessian, n, n);
- for (i = 1; i <= n; i++) {
+ for (i = 0; i < n; i++) {
tau_e[i]=x1[i] + x2[i];
}
- for (i = 1; i <= n; i++) {
+ for (i = 0; i < n; i++) {
tau_s[i]=x2[i];
}
- free_dvector(x1, 1, n);
- free_dvector(x2, 1, n);
+ free_dvector(x1, n);
+ free_dvector(x2, n);
}
void initialize(double f1, double f2, int n, double Q, double *x1, double *x2)
@@ -166,14 +166,14 @@ int i;
double q, omega, *tau_e, *tau_s;
double exp1, exp2, dexp, expo;
-tau_e = dvector(1, n);
-tau_s = dvector(1, n);
+tau_e = dvector(n);
+tau_s = dvector(n);
if (n > 1) {
exp1 = log10(f1);
exp2 = log10(f2);
dexp = (exp2 - exp1) / ((double) (n - 1));
q = 1.0 / ((n - 1.0) * Q);
- for (i = 1, expo = exp1; i <= n; i++, expo += dexp) {
+ for (i = 0, expo = exp1; i < n; i++, expo += dexp) {
omega = PI2 * pow(10.0, expo);
tau_s[i] = 1.0 / omega;
tau_e[i] = tau_s[i] * (1.0 + q) / (1.0 - q);
@@ -184,19 +184,19 @@ if (n > 1) {
exp2 = log10(f2);
expo=(exp1+exp2)/2.0;
omega = PI2 * pow(10.0, expo);
- tau_s[1] = 1.0 / omega;
- tau_e[1] = tau_s[1] * (1.0 + q) / (1.0 - q);
+ tau_s[0] = 1.0 / omega;
+ tau_e[0] = tau_s[0] * (1.0 + q) / (1.0 - q);
}
/*
* x1 denotes the parameter tau_e - tau_s and x2 denotes the parameter tau_s
*/
-for (i = 1; i <= n; i++) {
+for (i = 0; i < n; i++) {
x1[i] = tau_e[i] - tau_s[i];
x2[i] = tau_s[i];
}
-free_dvector(tau_e, 1, n);
-free_dvector(tau_s, 1, n);
+free_dvector(tau_e, n);
+free_dvector(tau_s, n);
}
void derivatives(double f1, double f2, int n, double Q, double *x1, double *x2, double *gradient, double **hessian)
@@ -207,15 +207,15 @@ double f, df, omega;
double *dadp, *dbdp, *dqdp, d2qdp2;
double tau_e, tau_s, a, b, Q_omega;
-dadp = dvector(1, n);
-dbdp = dvector(1, n);
-dqdp = dvector(1, n);
+dadp = dvector(n);
+dbdp = dvector(n);
+dqdp = dvector(n);
exp1 = log10(f1);
exp2 = log10(f2);
dexp = (exp2 - exp1) / 100.0;
-for (i = 1; i <= n; i++) {
+for (i = 0; i < n; i++) {
gradient[i] = 0.0;
- for (j = 1; j <= i; j++) {
+ for (j = 0; j < i; j++) {
hessian[j][i] = 0.0;
hessian[j][i] = hessian[i][j];
}
@@ -226,7 +226,7 @@ for (expo = exp1; expo <= exp2; expo += dexp) {
omega = PI2 * f;
a = (double) (1 - n);
b = 0.0;
- for (i = 1; i <= n; i++) {
+ for (i = 0; i < n; i++) {
tau_e = x1[i] + x2[i];
tau_s = x2[i];
a += (1.0 + omega * omega * tau_e * tau_s) /
@@ -237,10 +237,10 @@ for (expo = exp1; expo <= exp2; expo += dexp) {
dbdp[i] = omega / (1.0 + omega * omega * tau_s * tau_s);
}
Q_omega = a / b;
- for (i = 1; i <= n; i++) {
+ for (i = 0; i < n; i++) {
dqdp[i] = (dbdp[i] - (b / a) * dadp[i]) / a;
gradient[i] += 2.0 * (1.0 / Q_omega - 1.0 / Q) * dqdp[i] * df / (f2 - f1);
- for (j = 1; j <= i; j++) {
+ for (j = 0; j <= i; j++) {
d2qdp2 = -(dadp[i] * dbdp[j] + dbdp[i] * dadp[j]
- 2.0 * (b / a) * dadp[i] * dadp[j]) / (a * a);
hessian[i][j] += (2.0 * dqdp[i] * dqdp[j] + 2.0 * (1.0 / Q_omega - 1.0 / Q) * d2qdp2)
@@ -249,9 +249,9 @@ for (expo = exp1; expo <= exp2; expo += dexp) {
}
}
}
-free_dvector(dadp, 1, n);
-free_dvector(dbdp, 1, n);
-free_dvector(dqdp, 1, n);
+free_dvector(dadp, n);
+free_dvector(dbdp, n);
+free_dvector(dqdp, n);
}
void invert(double *x, double *b, double **A, int n)
@@ -259,32 +259,32 @@ void invert(double *x, double *b, double **A, int n)
int i, j, k;
double *xp, *W, **V, **A_inverse;
-xp = dvector(1, n);
-W = dvector(1, n);
-V = dmatrix(1, n, 1, n);
-A_inverse = dmatrix(1, n, 1, n);
+xp = dvector(n);
+W = dvector(n);
+V = dmatrix(n, n);
+A_inverse = dmatrix(n, n);
dsvdcmp(A, n, n, W, V);
-for (i = 1; i <= n; i++)
- for (j = 1; j <= n; j++)
+for (i = 0; i < n; i++)
+ for (j = 0; j < n; j++)
V[i][j] = (1.0 / W[i]) * A[j][i];
-for (i = 1; i <= n; i++) {
- for (j = 1; j <= n; j++) {
+for (i = 0; i < n; i++) {
+ for (j = 0; j < n; j++) {
A_inverse[i][j] = 0.0;
- for (k = 1; k <= n; k++)
+ for (k = 0; k < n; k++)
A_inverse[i][j] += A[i][k] * V[k][j];
}
}
-free_dvector(W, 1, n);
-free_dmatrix(V, 1, n, 1, n);
-for (i = 1; i <= n; i++) {
+free_dvector(W, n);
+free_dmatrix(V, n, n);
+for (i = 0; i < n; i++) {
xp[i] = x[i];
- for (j = 1; j <= n; j++) {
+ for (j = 0; j < n; j++) {
xp[i] -= A_inverse[i][j] * b[j];
}
x[i] = xp[i];
}
-free_dvector(xp, 1, n);
-free_dmatrix(A_inverse, 1, n, 1, n);
+free_dvector(xp, n);
+free_dmatrix(A_inverse, n, n);
}
void nrerror(const char *error_text)
@@ -295,43 +295,41 @@ void nrerror(const char *error_text)
exit(1);
}
-double *dvector(int nl, int nh)
+double *dvector(int n)
{
double *v;
- v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
+ v=(double *)malloc(n*sizeof(double));
if (!v) nrerror("allocation failure in dvector()");
- return v-nl;
+ return v;
}
-double **dmatrix(int nrl, int nrh, int ncl, int nch)
+double **dmatrix(int nr, int nc)
{
int i;
double **m;
- m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
+ m=(double **) malloc(nr*sizeof(double*));
if (!m) nrerror("allocation failure 1 in dmatrix()");
- m -= nrl;
- for(i=nrl;i<=nrh;i++) {
- m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
+ for(i=0;i<nr;i++) {
+ m[i]=(double *) malloc(nc*sizeof(double));
if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
- m[i] -= ncl;
}
return m;
}
-void free_dvector(double *v, int nl, int nh)
+void free_dvector(double *v, int n)
{
- free((char*) (v+nl));
+ free(v);
}
-void free_dmatrix(double **m, int nrl, int nrh, int ncl, int nch)
+void free_dmatrix(double **m, int nr, int nc)
{
int i;
- for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
- free((char*) (m+nrl));
+ for(i=0;i<nr;i++) free(m[i]);
+ free(m);
}
void dsvdcmp(double **a, int m, int n, double *w, double **v)
@@ -342,15 +340,15 @@ void dsvdcmp(double **a, int m, int n, double *w, double **v)
double *rv1;
if (m < n) nrerror("SVDCMP: You must augment A with extra zero rows");
- rv1=dvector(1,n);
- for (i=1;i<=n;i++) {
+ rv1=dvector(n);
+ for (i=0;i<n;i++) {
l=i+1;
rv1[i]=scale*g;
g=s=scale=0.0;
- if (i <= m) {
- for (k=i;k<=m;k++) scale += fabs(a[k][i]);
+ if (i < m) {
+ for (k=i;k<m;k++) scale += fabs(a[k][i]);
if (scale) {
- for (k=i;k<=m;k++) {
+ for (k=i;k<m;k++) {
a[k][i] /= scale;
s += a[k][i]*a[k][i];
}
@@ -358,22 +356,22 @@ void dsvdcmp(double **a, int m, int n, double *w, double **v)
g = -copysign(sqrt(s),f);
h=f*g-s;
a[i][i]=f-g;
- if (i != n) {
- for (j=l;j<=n;j++) {
- for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
+ if (i != n-1) {
+ for (j=l;j<n;j++) {
+ for (s=0.0,k=i;k<m;k++) s += a[k][i]*a[k][j];
f=s/h;
- for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
+ for (k=i;k<m;k++) a[k][j] += f*a[k][i];
}
}
- for (k=i;k<=m;k++) a[k][i] *= scale;
+ for (k=i;k<m;k++) a[k][i] *= scale;
}
}
w[i]=scale*g;
g=s=scale=0.0;
- if (i <= m && i != n) {
- for (k=l;k<=n;k++) scale += fabs(a[i][k]);
+ if (i < m && i != (n-1)) {
+ for (k=l;k<n;k++) scale += fabs(a[i][k]);
if (scale) {
- for (k=l;k<=n;k++) {
+ for (k=l;k<n;k++) {
a[i][k] /= scale;
s += a[i][k]*a[i][k];
}
@@ -381,58 +379,58 @@ void dsvdcmp(double **a, int m, int n, double *w, double **v)
g = -copysign(sqrt(s),f);
h=f*g-s;
a[i][l]=f-g;
- for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
- if (i != m) {
- for (j=l;j<=m;j++) {
- for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
- for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
+ for (k=l;k<n;k++) rv1[k]=a[i][k]/h;
+ if (i != (m-1)) {
+ for (j=l;j<m;j++) {
+ for (s=0.0,k=l;k<n;k++) s += a[j][k]*a[i][k];
+ for (k=l;k<n;k++) a[j][k] += s*rv1[k];
}
}
- for (k=l;k<=n;k++) a[i][k] *= scale;
+ for (k=l;k<n;k++) a[i][k] *= scale;
}
}
anorm=fmax(anorm, fabs(w[i])+fabs(rv1[i]));
}
- for (i=n;i>=1;i--) {
- if (i < n) {
+ for (i=n-1;i>=0;i--) {
+ if (i < n-1) {
if (g) {
- for (j=l;j<=n;j++)
+ for (j=l;j<n;j++)
v[j][i]=(a[i][j]/a[i][l])/g;
- for (j=l;j<=n;j++) {
- for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
- for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
+ for (j=l;j<n;j++) {
+ for (s=0.0,k=l;k<n;k++) s += a[i][k]*v[k][j];
+ for (k=l;k<n;k++) v[k][j] += s*v[k][i];
}
}
- for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
+ for (j=l;j<n;j++) v[i][j]=v[j][i]=0.0;
}
v[i][i]=1.0;
g=rv1[i];
l=i;
}
- for (i=n;i>=1;i--) {
+ for (i=n-1;i>=0;i--) {
l=i+1;
g=w[i];
- if (i < n)
- for (j=l;j<=n;j++) a[i][j]=0.0;
+ if (i < n-1)
+ for (j=l;j<n;j++) a[i][j]=0.0;
if (g) {
g=1.0/g;
- if (i != n) {
- for (j=l;j<=n;j++) {
- for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
+ if (i != n-1) {
+ for (j=l;j<n;j++) {
+ for (s=0.0,k=l;k<m;k++) s += a[k][i]*a[k][j];
f=(s/a[i][i])*g;
- for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
+ for (k=i;k<m;k++) a[k][j] += f*a[k][i];
}
}
- for (j=i;j<=m;j++) a[j][i] *= g;
+ for (j=i;j<m;j++) a[j][i] *= g;
} else {
- for (j=i;j<=m;j++) a[j][i]=0.0;
+ for (j=i;j<m;j++) a[j][i]=0.0;
}
++a[i][i];
}
- for (k=n;k>=1;k--) {
- for (its=1;its<=30;its++) {
+ for (k=n-1;k>=0;k--) {
+ for (its=0;its<30;its++) {
flag=1;
- for (l=k;l>=1;l--) {
+ for (l=k;l>=0;l--) {
nm=l-1;
if (fabs(rv1[l])+anorm == anorm) {
flag=0;
@@ -443,7 +441,7 @@ void dsvdcmp(double **a, int m, int n, double *w, double **v)
if (flag) {
c=0.0;
s=1.0;
- for (i=l;i<=k;i++) {
+ for (i=l;i<k;i++) {
f=s*rv1[i];
if (fabs(f)+anorm != anorm) {
g=w[i];
@@ -452,7 +450,7 @@ void dsvdcmp(double **a, int m, int n, double *w, double **v)
h=1.0/h;
c=g*h;
s=(-f*h);
- for (j=1;j<=m;j++) {
+ for (j=0;j<m;j++) {
y=a[j][nm];
z=a[j][i];
a[j][nm]=y*c+z*s;
@@ -465,7 +463,7 @@ void dsvdcmp(double **a, int m, int n, double *w, double **v)
if (l == k) {
if (z < 0.0) {
w[k] = -z;
- for (j=1;j<=n;j++) v[j][k]=(-v[j][k]);
+ for (j=0;j<n;j++) v[j][k]=(-v[j][k]);
}
break;
}
@@ -493,7 +491,7 @@ void dsvdcmp(double **a, int m, int n, double *w, double **v)
g=g*c-x*s;
h=y*s;
y=y*c;
- for (jj=1;jj<=n;jj++) {
+ for (jj=0;jj<n;jj++) {
x=v[jj][j];
z=v[jj][i];
v[jj][j]=x*c+z*s;
@@ -508,7 +506,7 @@ void dsvdcmp(double **a, int m, int n, double *w, double **v)
}
f=(c*g)+(s*y);
x=(c*y)-(s*g);
- for (jj=1;jj<=m;jj++) {
+ for (jj=0;jj<m;jj++) {
y=a[jj][j];
z=a[jj][i];
a[jj][j]=y*c+z*s;
@@ -520,6 +518,6 @@ void dsvdcmp(double **a, int m, int n, double *w, double **v)
w[k]=x;
}
}
- free_dvector(rv1,1,n);
+ free_dvector(rv1,n);
}
More information about the CIG-COMMITS
mailing list