Annotation of OpenXM/src/asir-contrib/packages/doc/nk_fb_gen_c/test1.c, Revision 1.1
1.1 ! takayama 1: /* $OpenXM$ */
! 2: #include "test1.h"
! 3: #define SEARCH_MIN_ITER 100
! 4: int main()
! 5: {
! 6: extern double g_y1;
! 7: extern double g_y2;
! 8: extern double g_x11;
! 9: extern double g_x12;
! 10: extern double g_x22;
! 11: extern double g_s1;
! 12: extern double g_s2;
! 13: extern double g_r;
! 14: extern double g_s11;
! 15: extern double g_s12;
! 16: extern double g_s22;
! 17: extern double values[N_VALUES];
! 18:
! 19: init_mat();
! 20: g_r = 1.0;
! 21: /* Write data here. */
! 22: /* g_x11, g_x12, g_x22, g_y1, g_y2 */
! 23: /* g_s11, g_s12, g_s22, g_s1, g_s2 */
! 24: /* Wind data from this11/h-mle/FB2/s1_wind.c
! 25: The data average s1 and s2 is defined by
! 26: (s1,s2) = (sum(s1_i)/N,sum(s2_i)/N)
! 27: where (s1_i,s2_i), i=1,2, ..., N are data on the S^1 (circle).
! 28: s11 is defined as the average of s1_i*s1_i; s11=sum(s1_i*s1_i)/N.
! 29: s12 is defined as the average of s1_i*s2_i; s11=sum(s1_i*s2_i)/N.
! 30: s22 is defined as the average of s2_i*s2_i; s11=sum(s2_i*s2_i)/N.
! 31:
! 32: */
! 33: g_s11 = 0.662125; g_s12 = 0.274563; g_s22 = 0.337875; g_s1 = 0.317564; g_s2 = -0.020188;
! 34: /* Initial parameter for the HGD. */
! 35: g_x11 = 0.5; g_x12 = 0.5; g_x22 = 0.15; g_y1 = 1; g_y2 = 0.15; g_r = 1.0;
! 36:
! 37: double *val;
! 38: val = init_val(1, g_x11, g_x12, g_x22, g_y1, g_y2, g_s11, g_s12, g_s22, g_s1, g_s2);
! 39: show_v(val, N_VALUES);
! 40: search_min(val);
! 41: printf("search_min :\n");
! 42: show_v(values, N_VALUES);
! 43: val = init_val(1, g_x11, g_x12, g_x22, g_y1, g_y2, g_s11, g_s12, g_s22, g_s1, g_s2);
! 44: printf("init_val :\n");
! 45: show_v(val, N_VALUES);
! 46: }
! 47:
! 48: void init_mat()
! 49: {
! 50: extern gsl_matrix *a1, *b1, *c1, *e1;
! 51: extern gsl_matrix *a2, *b2, *c2, *e2;
! 52: extern gsl_matrix *a3, *b3, *c3, *e3;
! 53: extern gsl_matrix *p2, *q2, *p3, *q3, *r3;
! 54: extern gsl_matrix *db11;
! 55: extern gsl_matrix *dc11;
! 56: extern gsl_matrix *db12;
! 57: extern gsl_matrix *dc12;
! 58: extern gsl_matrix *db21;
! 59: extern gsl_matrix *dc21;
! 60: extern gsl_matrix *db22;
! 61: extern gsl_matrix *dc22;
! 62: extern gsl_matrix *dq21;
! 63: extern gsl_matrix *dq31;
! 64: extern gsl_matrix *dr31;
! 65: extern gsl_matrix *dq22;
! 66: extern gsl_matrix *dq32;
! 67: extern gsl_matrix *dr32;
! 68: extern gsl_matrix *inv_a1;
! 69: extern gsl_matrix *inv_a2;
! 70: extern gsl_matrix *inv_a3;
! 71: extern gsl_matrix *inv_p2;
! 72: extern gsl_matrix *inv_p3;
! 73: extern gsl_matrix *pf1_m;
! 74: extern gsl_matrix *pf2_m;
! 75: extern gsl_matrix *pf1_nd_m;
! 76: extern gsl_matrix *pf2_nd_m;
! 77: extern gsl_matrix *pf11_m;
! 78: extern gsl_matrix *pf12_m;
! 79: extern gsl_matrix *pf22_m;
! 80: extern gsl_matrix *pft_m;
! 81: extern gsl_vector *grad_v;
! 82:
! 83: a1 = gsl_matrix_alloc(4, 4);
! 84: b1 = gsl_matrix_alloc(4, 4);
! 85: c1 = gsl_matrix_alloc(4, 1);
! 86: e1 = gsl_matrix_alloc(4, 2);
! 87: a2 = gsl_matrix_alloc(4, 4);
! 88: b2 = gsl_matrix_alloc(4, 4);
! 89: c2 = gsl_matrix_alloc(4, 1);
! 90: e2 = gsl_matrix_alloc(4, 2);
! 91: p2 = gsl_matrix_alloc(1, 1);
! 92: q2 = gsl_matrix_alloc(1, 4);
! 93: p3 = gsl_matrix_alloc(2, 2);
! 94: q3 = gsl_matrix_alloc(2, 1);
! 95: r3 = gsl_matrix_alloc(2, 4);
! 96: db11 = gsl_matrix_alloc(4, 4);
! 97: dc11 = gsl_matrix_alloc(4, 1);
! 98: db12 = gsl_matrix_alloc(4, 4);
! 99: dc12 = gsl_matrix_alloc(4, 1);
! 100: db21 = gsl_matrix_alloc(4, 4);
! 101: dc21 = gsl_matrix_alloc(4, 1);
! 102: db22 = gsl_matrix_alloc(4, 4);
! 103: dc22 = gsl_matrix_alloc(4, 1);
! 104: dq21 = gsl_matrix_alloc(1, 4);
! 105: dq31 = gsl_matrix_alloc(2, 1);
! 106: dr31 = gsl_matrix_alloc(2, 4);
! 107: dq22 = gsl_matrix_alloc(1, 4);
! 108: dq32 = gsl_matrix_alloc(2, 1);
! 109: dr32 = gsl_matrix_alloc(2, 4);
! 110: inv_a1 = gsl_matrix_alloc(4, 4);
! 111: inv_a2 = gsl_matrix_alloc(4, 4);
! 112: inv_p2 = gsl_matrix_alloc(1, 1);
! 113: inv_p3 = gsl_matrix_alloc(2, 2);
! 114: pf1_m = gsl_matrix_alloc(4, 4);
! 115: pf2_m = gsl_matrix_alloc(4, 4);
! 116: pf1_nd_m = gsl_matrix_alloc(4, 4);
! 117: pf2_nd_m = gsl_matrix_alloc(4, 4);
! 118: pf11_m = gsl_matrix_alloc(4, 4);
! 119: pf12_m = gsl_matrix_alloc(4, 4);
! 120: pf22_m = gsl_matrix_alloc(4, 4);
! 121: pft_m = gsl_matrix_alloc(N_VALUES, N_VALUES);
! 122: grad_v = gsl_vector_alloc(DIM);
! 123: }
! 124:
! 125: void set_abce_1(double x11, double x12, double x22, double y1, double y2, double r)
! 126: {
! 127: extern gsl_matrix *a1, *b1, *c1, *e1;
! 128: gsl_matrix_set(a1, 0, 0, 1);
! 129: gsl_matrix_set(a1, 1, 1, 1);
! 130: gsl_matrix_set(a1, 2, 1, 1*x12);
! 131: gsl_matrix_set(a1, 2, 2, (-2)*x11 + 2*x22);
! 132: gsl_matrix_set(a1, 3, 3, 1);
! 133: gsl_matrix_set(b1, 0, 1, 1);
! 134: gsl_matrix_set(b1, 1, 3, 1);
! 135: gsl_matrix_set(b1, 2, 0, 1*x12*r*r);
! 136: gsl_matrix_set(b1, 2, 1, (-1)*y2);
! 137: gsl_matrix_set(b1, 2, 2, 1*y1);
! 138: gsl_matrix_set(b1, 2, 3, (-1)*x12);
! 139: gsl_matrix_set(e1, 3, 0, 1);
! 140: }
! 141:
! 142: void set_abce_2(double x11, double x12, double x22, double y1, double y2, double r)
! 143: {
! 144: extern gsl_matrix *a2, *b2, *c2, *e2;
! 145: gsl_matrix_set(a2, 0, 0, 1);
! 146: gsl_matrix_set(a2, 1, 1, 2*x11 + (-2)*x22);
! 147: gsl_matrix_set(a2, 1, 2, 1*x12);
! 148: gsl_matrix_set(a2, 2, 2, 1);
! 149: gsl_matrix_set(a2, 3, 3, (-2)*x11 + 2*x22);
! 150: gsl_matrix_set(b2, 0, 2, 1);
! 151: gsl_matrix_set(b2, 1, 1, 1*y2);
! 152: gsl_matrix_set(b2, 1, 2, (-1)*y1);
! 153: gsl_matrix_set(b2, 1, 3, 1*x12);
! 154: gsl_matrix_set(b2, 2, 0, 1*r*r);
! 155: gsl_matrix_set(b2, 2, 3, (-1));
! 156: gsl_matrix_set(b2, 3, 1, 1*x12*r*r);
! 157: gsl_matrix_set(b2, 3, 2, 1);
! 158: gsl_matrix_set(b2, 3, 3, (-1)*y2);
! 159: gsl_matrix_set(c2, 3, 0, 1*y1);
! 160: gsl_matrix_set(e2, 3, 0, (-2)*x12);
! 161: }
! 162:
! 163: void set_pqr(double x11, double x12, double x22, double y1, double y2, double r)
! 164: {
! 165: extern gsl_matrix *p2, *q2, *p3, *q3, *r3;
! 166: gsl_matrix_set(p2, 0, 0, (-2)*x11 + 2*x22);
! 167: gsl_matrix_set(q2, 0, 0, (-1)*x12*r*r);
! 168: gsl_matrix_set(q2, 0, 1, 1*y2);
! 169: gsl_matrix_set(q2, 0, 2, (-1)*y1);
! 170: gsl_matrix_set(q2, 0, 3, 2*x12);
! 171: gsl_matrix_set(p3, 0, 0, 2*x11 + (-2)*x22);
! 172: gsl_matrix_set(p3, 0, 1, 2*x12);
! 173: gsl_matrix_set(p3, 1, 0, 2*x12);
! 174: gsl_matrix_set(p3, 1, 1, (-2)*x11 + 2*x22);
! 175: gsl_matrix_set(q3, 0, 0, 1*y2);
! 176: gsl_matrix_set(q3, 1, 0, (-1)*y1);
! 177: gsl_matrix_set(r3, 0, 0, (-1)*y1*r*r);
! 178: gsl_matrix_set(r3, 0, 1, (-2)*x11*r*r + 2*r*r*x22 + 1);
! 179: gsl_matrix_set(r3, 0, 2, (-1)*x12*r*r);
! 180: gsl_matrix_set(r3, 0, 3, 1*y1);
! 181: gsl_matrix_set(r3, 1, 1, (-1)*x12*r*r);
! 182: gsl_matrix_set(r3, 1, 2, (-1));
! 183: gsl_matrix_set(r3, 1, 3, 1*y2);
! 184: }
! 185:
! 186: void set_dbcqr(double x11, double x12, double x22, double y1, double y2, double r)
! 187: {
! 188: extern gsl_matrix *db11;
! 189: extern gsl_matrix *dc11;
! 190: extern gsl_matrix *db12;
! 191: extern gsl_matrix *dc12;
! 192: extern gsl_matrix *db21;
! 193: extern gsl_matrix *dc21;
! 194: extern gsl_matrix *db22;
! 195: extern gsl_matrix *dc22;
! 196: extern gsl_matrix *dq21;
! 197: extern gsl_matrix *dq31;
! 198: extern gsl_matrix *dr31;
! 199: extern gsl_matrix *dq22;
! 200: extern gsl_matrix *dq32;
! 201: extern gsl_matrix *dr32;
! 202: gsl_matrix_set(db11, 2, 2, 1);
! 203: gsl_matrix_set(db12, 2, 1, (-1));
! 204: gsl_matrix_set(db21, 1, 2, (-1));
! 205: gsl_matrix_set(dc21, 3, 0, 1);
! 206: gsl_matrix_set(db22, 1, 1, 1);
! 207: gsl_matrix_set(db22, 3, 3, (-1));
! 208: gsl_matrix_set(dq21, 0, 2, (-1));
! 209: gsl_matrix_set(dq31, 1, 0, (-1));
! 210: gsl_matrix_set(dr31, 0, 0, (-1)*r*r);
! 211: gsl_matrix_set(dr31, 0, 3, 1);
! 212: gsl_matrix_set(dq22, 0, 1, 1);
! 213: gsl_matrix_set(dq32, 0, 0, 1);
! 214: gsl_matrix_set(dr32, 1, 3, 1);
! 215: }
! 216:
! 217: void pf_all(double x11, double x12, double x22, double y1, double y2, double r)
! 218: {
! 219: extern gsl_matrix *a1, *inv_a1, *b1, *c1, *e1;
! 220: extern gsl_matrix *a2, *inv_a2, *b2, *c2, *e2;
! 221: extern gsl_matrix *p2, *inv_p2, *q2;
! 222: extern gsl_matrix *p3, *inv_p3, *q3, *r3;
! 223: extern gsl_matrix *pf1_m;
! 224: extern gsl_matrix *pf1_nd_m;
! 225: extern gsl_matrix *pf2_m;
! 226: extern gsl_matrix *pf2_nd_m;
! 227: extern double g_s1;
! 228: extern double g_s2;
! 229:
! 230: extern gsl_matrix *db11;
! 231: extern gsl_matrix *dc11;
! 232: extern gsl_matrix *db12;
! 233: extern gsl_matrix *dc12;
! 234: extern gsl_matrix *db21;
! 235: extern gsl_matrix *dc21;
! 236: extern gsl_matrix *db22;
! 237: extern gsl_matrix *dc22;
! 238: extern gsl_matrix *dq21;
! 239: extern gsl_matrix *dq31;
! 240: extern gsl_matrix *dr31;
! 241: extern gsl_matrix *dq22;
! 242: extern gsl_matrix *dq32;
! 243: extern gsl_matrix *dr32;
! 244:
! 245: extern gsl_matrix *pf11_m;
! 246: extern gsl_matrix *pf12_m;
! 247: extern gsl_matrix *pf22_m;
! 248: extern double g_s11;
! 249: extern double g_s12;
! 250: extern double g_s22;
! 251: gsl_matrix *t1;
! 252: gsl_matrix *t2;
! 253: gsl_matrix *t3;
! 254: gsl_matrix *t4;
! 255: gsl_matrix *t5;
! 256: gsl_matrix *t6;
! 257:
! 258: static double o_y1;
! 259: static double o_y2;
! 260: static double o_x11;
! 261: static double o_x12;
! 262: static double o_x22;
! 263: if (o_y1 == y1 && o_y2 == y2 && o_x11 == x11 &&o_x12 == x12 &&o_x22 == x22)
! 264: return;
! 265: o_y1 = y1;
! 266: o_y2 = y2;
! 267: o_x11 = x11;
! 268: o_x12 = x12;
! 269: o_x22 = x22;
! 270:
! 271: set_pqr(x11, x12, x22, y1, y2, r);
! 272: invmat(p2, inv_p2);
! 273: invmat(p3, inv_p3);
! 274:
! 275: set_abce_1(x11, x12, x22, y1, y2, r);
! 276: invmat(a1, inv_a1);
! 277: t1 = gsl_matrix_alloc(inv_p2->size1, q2->size2);
! 278: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p2, q2, 0.0, t1);
! 279: t2 = gsl_matrix_alloc(c1->size1, t1->size2);
! 280: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c1, t1, 0.0, t2);
! 281: t3 = gsl_matrix_alloc(q3->size1, t1->size2);
! 282: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, q3, t1, 0.0, t3);
! 283: gsl_matrix_sub(t3, r3);
! 284: t4 = gsl_matrix_alloc(inv_p3->size1, t3->size2);
! 285: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p3, t3, 0.0, t4);
! 286: t5 = gsl_matrix_alloc(e1->size1, t4->size2);
! 287: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, e1, t4, 0.0, t5);
! 288: gsl_matrix_add(t5, b1);
! 289: gsl_matrix_sub(t5, t2);
! 290: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_a1, t5, 0.0, t2);
! 291: gsl_matrix_memcpy(pf1_nd_m, t2);
! 292: gsl_matrix_memcpy(pf1_m, t2);
! 293:
! 294: gsl_matrix_set_identity(t2);
! 295: gsl_matrix_scale(t2, -g_s1);
! 296: gsl_matrix_add(pf1_m, t2);
! 297:
! 298: gsl_matrix_free(t1);
! 299: gsl_matrix_free(t2);
! 300: gsl_matrix_free(t3);
! 301: gsl_matrix_free(t4);
! 302: gsl_matrix_free(t5);
! 303:
! 304: set_abce_2(x11, x12, x22, y1, y2, r);
! 305: invmat(a2, inv_a2);
! 306: t1 = gsl_matrix_alloc(inv_p2->size1, q2->size2);
! 307: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p2, q2, 0.0, t1);
! 308: t2 = gsl_matrix_alloc(c2->size1, t1->size2);
! 309: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c2, t1, 0.0, t2);
! 310: t3 = gsl_matrix_alloc(q3->size1, t1->size2);
! 311: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, q3, t1, 0.0, t3);
! 312: gsl_matrix_sub(t3, r3);
! 313: t4 = gsl_matrix_alloc(inv_p3->size1, t3->size2);
! 314: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p3, t3, 0.0, t4);
! 315: t5 = gsl_matrix_alloc(e2->size1, t4->size2);
! 316: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, e2, t4, 0.0, t5);
! 317: gsl_matrix_add(t5, b2);
! 318: gsl_matrix_sub(t5, t2);
! 319: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_a2, t5, 0.0, t2);
! 320: gsl_matrix_memcpy(pf2_nd_m, t2);
! 321: gsl_matrix_memcpy(pf2_m, t2);
! 322:
! 323: gsl_matrix_set_identity(t2);
! 324: gsl_matrix_scale(t2, -g_s2);
! 325: gsl_matrix_add(pf2_m, t2);
! 326:
! 327: gsl_matrix_free(t1);
! 328: gsl_matrix_free(t2);
! 329: gsl_matrix_free(t3);
! 330: gsl_matrix_free(t4);
! 331: gsl_matrix_free(t5);
! 332:
! 333: set_dbcqr(x11, x12, x22, y1, y2, r);
! 334: t1 = gsl_matrix_alloc(inv_p2->size1, q2->size2);
! 335: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p2, q2, 0.0, t1);
! 336: t2 = gsl_matrix_alloc(inv_p2->size1, dq21->size2);
! 337: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p2, dq21, 0.0, t2);
! 338: t3 = gsl_matrix_alloc(dq31->size1, t1->size2);
! 339: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, dq31, t1, 0.0, t3);
! 340: t4 = gsl_matrix_alloc(q3->size1, t2->size2);
! 341: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, q3, t2, 0.0, t4);
! 342: gsl_matrix_add(t3, t4);
! 343: gsl_matrix_sub(t3, dr31);
! 344: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p3, t3, 0.0, t4);
! 345: t5 = gsl_matrix_alloc(e1->size1, t4->size2);
! 346: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, e1, t4, 0.0, t5);
! 347: gsl_matrix_add(t5, db11);
! 348: t6 = gsl_matrix_alloc(dc11->size1, t1->size2);
! 349: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, dc11, t1, 0.0, t6);
! 350: gsl_matrix_sub(t5, t6);
! 351: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c1, t2, 0.0, t6);
! 352: gsl_matrix_sub(t5, t6);
! 353: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_a1, t5, 0.0, t6);
! 354: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, pf1_nd_m, pf1_nd_m, 0.0, t5);
! 355: gsl_matrix_add(t5, t6);
! 356: gsl_matrix_memcpy(pf11_m, t5);
! 357:
! 358: gsl_matrix_set_identity(t5);
! 359: gsl_matrix_scale(t5, -g_s11);
! 360: gsl_matrix_add(pf11_m, t5);
! 361: gsl_matrix_free(t1);
! 362: gsl_matrix_free(t2);
! 363: gsl_matrix_free(t3);
! 364: gsl_matrix_free(t4);
! 365: gsl_matrix_free(t5);
! 366: gsl_matrix_free(t6);
! 367:
! 368: t1 = gsl_matrix_alloc(inv_p2->size1, q2->size2);
! 369: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p2, q2, 0.0, t1);
! 370: t2 = gsl_matrix_alloc(inv_p2->size1, dq22->size2);
! 371: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p2, dq22, 0.0, t2);
! 372: t3 = gsl_matrix_alloc(dq32->size1, t1->size2);
! 373: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, dq32, t1, 0.0, t3);
! 374: t4 = gsl_matrix_alloc(q3->size1, t2->size2);
! 375: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, q3, t2, 0.0, t4);
! 376: gsl_matrix_add(t3, t4);
! 377: gsl_matrix_sub(t3, dr32);
! 378: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p3, t3, 0.0, t4);
! 379: t5 = gsl_matrix_alloc(e1->size1, t4->size2);
! 380: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, e1, t4, 0.0, t5);
! 381: gsl_matrix_add(t5, db12);
! 382: t6 = gsl_matrix_alloc(dc12->size1, t1->size2);
! 383: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, dc12, t1, 0.0, t6);
! 384: gsl_matrix_sub(t5, t6);
! 385: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c1, t2, 0.0, t6);
! 386: gsl_matrix_sub(t5, t6);
! 387: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_a1, t5, 0.0, t6);
! 388: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, pf1_nd_m, pf2_nd_m, 0.0, t5);
! 389: gsl_matrix_add(t5, t6);
! 390: gsl_matrix_memcpy(pf12_m, t5);
! 391:
! 392: gsl_matrix_set_identity(t5);
! 393: gsl_matrix_scale(t5, -g_s12);
! 394: gsl_matrix_add(pf12_m, t5);
! 395: gsl_matrix_free(t1);
! 396: gsl_matrix_free(t2);
! 397: gsl_matrix_free(t3);
! 398: gsl_matrix_free(t4);
! 399: gsl_matrix_free(t5);
! 400: gsl_matrix_free(t6);
! 401:
! 402: t1 = gsl_matrix_alloc(inv_p2->size1, q2->size2);
! 403: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p2, q2, 0.0, t1);
! 404: t2 = gsl_matrix_alloc(inv_p2->size1, dq22->size2);
! 405: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p2, dq22, 0.0, t2);
! 406: t3 = gsl_matrix_alloc(dq32->size1, t1->size2);
! 407: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, dq32, t1, 0.0, t3);
! 408: t4 = gsl_matrix_alloc(q3->size1, t2->size2);
! 409: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, q3, t2, 0.0, t4);
! 410: gsl_matrix_add(t3, t4);
! 411: gsl_matrix_sub(t3, dr32);
! 412: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_p3, t3, 0.0, t4);
! 413: t5 = gsl_matrix_alloc(e2->size1, t4->size2);
! 414: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, e2, t4, 0.0, t5);
! 415: gsl_matrix_add(t5, db22);
! 416: t6 = gsl_matrix_alloc(dc22->size1, t1->size2);
! 417: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, dc22, t1, 0.0, t6);
! 418: gsl_matrix_sub(t5, t6);
! 419: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c2, t2, 0.0, t6);
! 420: gsl_matrix_sub(t5, t6);
! 421: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, inv_a2, t5, 0.0, t6);
! 422: gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, pf2_nd_m, pf2_nd_m, 0.0, t5);
! 423: gsl_matrix_add(t5, t6);
! 424: gsl_matrix_memcpy(pf22_m, t5);
! 425:
! 426: gsl_matrix_set_identity(t5);
! 427: gsl_matrix_scale(t5, -g_s22);
! 428: gsl_matrix_add(pf22_m, t5);
! 429: gsl_matrix_free(t1);
! 430: gsl_matrix_free(t2);
! 431: gsl_matrix_free(t3);
! 432: gsl_matrix_free(t4);
! 433: gsl_matrix_free(t5);
! 434: gsl_matrix_free(t6);
! 435:
! 436: }
! 437:
! 438: void invmat(gsl_matrix *m, gsl_matrix *invm)
! 439: {
! 440: int n = m->size1;
! 441: int s;
! 442: gsl_permutation *p = gsl_permutation_alloc(n);
! 443: gsl_matrix *old_m = gsl_matrix_alloc(m->size1, m->size2);
! 444: gsl_matrix_memcpy(old_m, m);
! 445:
! 446: gsl_linalg_LU_decomp(m, p, &s);
! 447: gsl_linalg_LU_invert(m, p, invm);
! 448:
! 449: gsl_permutation_free(p);
! 450: gsl_matrix_memcpy(m, old_m);
! 451: gsl_matrix_free(old_m);
! 452: }
! 453:
! 454: int move_t(double x11, double x12, double x22, double y1, double y2, double xx11, double xx12, double xx22, double yy1, double yy2, double *val)
! 455: {
! 456: extern double g_x11;
! 457: extern double g_x12;
! 458: extern double g_x22;
! 459: extern double g_y1;
! 460: extern double g_y2;
! 461: extern double move_t_points[DIM];
! 462:
! 463: double params[DIM];
! 464: params[0] = xx11 - x11;
! 465: params[1] = xx12 - x12;
! 466: params[2] = xx22 - x22;
! 467: params[3] = yy1 - y1;
! 468: params[4] = yy2 - y2;
! 469:
! 470: const gsl_odeiv_step_type *T = ODEIV_STEP_TYPE;
! 471: gsl_odeiv_step *s = gsl_odeiv_step_alloc(T, N_VALUES);
! 472: gsl_odeiv_control *c = gsl_odeiv_control_y_new(1e-6, 0.0);
! 473: gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc(N_VALUES);
! 474: gsl_odeiv_system sys = {sys_t, NULL, N_VALUES, params};
! 475:
! 476: double t = 0.0;
! 477: double h = 1e-6;
! 478: double hmin = 1e-6;
! 479:
! 480: #ifdef DEBUG_RK
! 481: printf("Runge-Kutta: [%f,%f,%f,%f,%f] -> [%f,%f,%f,%f,%f]\n",x11, x12, x22, y1, y2, xx11, xx12, xx22, yy1, yy2);
! 482: #endif
! 483: while (t < 1.0) {
! 484: int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &t, 1.0, &h, val);
! 485: if (status != GSL_SUCCESS) break;
! 486: if (h < hmin) {
! 487: printf("move_t : h is too small\n");
! 488: move_t_points[0] = g_x11 + t * params[0];
! 489: move_t_points[1] = g_x12 + t * params[1];
! 490: move_t_points[2] = g_x22 + t * params[2];
! 491: move_t_points[3] = g_y1 + t * params[3];
! 492: move_t_points[4] = g_y2 + t * params[4];
! 493: gsl_odeiv_evolve_free(e);
! 494: gsl_odeiv_control_free(c);
! 495: gsl_odeiv_step_free(s);
! 496: return MOVE_T_FAIL;
! 497: }
! 498: #ifdef DEBUG_RK
! 499: printf("t = %f [%f,%f,%f,%f] : ", t, val[0],val[1],val[2],val[3]);
! 500: printf("[%f,%f,%f,%f,%f]\n", g_x11 + t * params[0], g_x12 + t * params[1], g_x22 + t * params[2], g_y1 + t * params[3], g_y2 + t * params[4]);
! 501: #endif
! 502: }
! 503: gsl_odeiv_evolve_free(e);
! 504: gsl_odeiv_control_free(c);
! 505: gsl_odeiv_step_free(s);
! 506: return MOVE_T_SUCCESS;
! 507: }
! 508:
! 509: int sys_t(double t, const double *y, double *val, double *params)
! 510: {
! 511: extern double g_x11;
! 512: extern double g_x12;
! 513: extern double g_x22;
! 514: extern double g_y1;
! 515: extern double g_y2;
! 516: extern double g_r;
! 517: extern gsl_matrix *pf11_m;
! 518: extern gsl_matrix *pf12_m;
! 519: extern gsl_matrix *pf22_m;
! 520: extern gsl_matrix *pf1_m;
! 521: extern gsl_matrix *pf2_m;
! 522: extern gsl_matrix *pft_m;
! 523:
! 524: int i, j;
! 525: double n_x11, n_x12, n_x22, n_y1, n_y2, n_r;
! 526:
! 527: gsl_matrix *t_m = gsl_matrix_alloc(pf11_m->size1, pf11_m->size2);
! 528:
! 529: n_x11 = g_x11 + t * params[0];
! 530: n_x12 = g_x12 + t * params[1];
! 531: n_x22 = g_x22 + t * params[2];
! 532: n_y1 = g_y1 + t * params[3];
! 533: n_y2 = g_y2 + t * params[4];
! 534: n_r = g_r;
! 535: gsl_matrix_set_zero(pft_m);
! 536: pf_all(n_x11, n_x12, n_x22, n_y1, n_y2, n_r);
! 537: gsl_matrix_memcpy(t_m, pf11_m);
! 538: gsl_matrix_scale(t_m, params[0]);
! 539: gsl_matrix_add(pft_m, t_m);
! 540: gsl_matrix_memcpy(t_m, pf12_m);
! 541: gsl_matrix_scale(t_m, params[1]);
! 542: gsl_matrix_add(pft_m, t_m);
! 543: gsl_matrix_memcpy(t_m, pf22_m);
! 544: gsl_matrix_scale(t_m, params[2]);
! 545: gsl_matrix_add(pft_m, t_m);
! 546: gsl_matrix_memcpy(t_m, pf1_m);
! 547: gsl_matrix_scale(t_m, params[3]);
! 548: gsl_matrix_add(pft_m, t_m);
! 549: gsl_matrix_memcpy(t_m, pf2_m);
! 550: gsl_matrix_scale(t_m, params[4]);
! 551: gsl_matrix_add(pft_m, t_m);
! 552: for (i = 0; i < N_VALUES; i++) {
! 553: val[i] = 0.0;
! 554: for (j = 0; j < N_VALUES; j++)
! 555: val[i] += gsl_matrix_get(pft_m, i, j) * y[j];
! 556: }
! 557: gsl_matrix_free(t_m);
! 558: return GSL_SUCCESS;
! 559: }
! 560:
! 561: double *init_val(int dim, double x11, double x12, double x22, double y1, double y2, double s11, double s12, double s22, double s1, double s2)
! 562: {
! 563: double xxx[MAXSIZE][MAXSIZE];
! 564: xxx[0][0] = x11;
! 565: xxx[0][1] = x12/2;
! 566: xxx[1][0] = x12/2;
! 567: xxx[1][1] = x22;
! 568: double y[] = {y1, y2};
! 569: double exp_part = exp(-(s11 * x11 + s12 * x12 + s22 * x22 + s1 * y1 + s2 * y2));
! 570: int maxdeg = 10;
! 571: int weight[] = {1, 1, 1, 1};
! 572: double *val = fbnd(dim, xxx, y, maxdeg, weight);
! 573: int i;
! 574: for(i = 0; i < 2*dim+2; i++)
! 575: val[i] = exp_part * val[i];
! 576: return val;
! 577: }
! 578:
! 579: void search_min(double *val)
! 580: {
! 581: extern double g_x11;
! 582: extern double g_x12;
! 583: extern double g_x22;
! 584: extern double g_y1;
! 585: extern double g_y2;
! 586: extern double values[N_VALUES];
! 587: extern gsl_vector *grad_v;
! 588:
! 589: size_t iter = 0;
! 590: int status;
! 591: const gsl_multimin_fdfminimizer_type *T;
! 592: gsl_multimin_fdfminimizer *s;
! 593: int move_t_status;
! 594: double t_x11;
! 595: double t_x12;
! 596: double t_x22;
! 597: double t_y1;
! 598: double t_y2;
! 599:
! 600: int j;
! 601: for (j = 0; j < N_VALUES; j++)
! 602: values[j] = val[j];
! 603: gsl_vector *x;
! 604: gsl_multimin_function_fdf my_func;
! 605: my_func.n = 5;
! 606: my_func.f = &my_f;
! 607: my_func.df = &my_df;
! 608: my_func.fdf = &my_fdf;
! 609: my_func.params = NULL;
! 610: x = gsl_vector_alloc(5);
! 611: gsl_vector_set(x, 0, g_x11);
! 612: gsl_vector_set(x, 1, g_x12);
! 613: gsl_vector_set(x, 2, g_x22);
! 614: gsl_vector_set(x, 3, g_y1);
! 615: gsl_vector_set(x, 4, g_y2);
! 616: T = MULTIMIN_FDFMINIMIZER_TYPE;
! 617: s = gsl_multimin_fdfminimizer_alloc(T, DIM);
! 618: gsl_multimin_fdfminimizer_set(s, &my_func, x, 0.01, 1e-4);
! 619: do {
! 620: iter++;
! 621: status = gsl_multimin_fdfminimizer_iterate(s);
! 622: if (status) {
! 623: if (status == GSL_ENOPROG)
! 624: printf("GSL_ENOPROG : gsl_multimin_fdminimizer_iterate\n");
! 625: break;
! 626: }
! 627: status = gsl_multimin_test_gradient(s->gradient, 1e-3);
! 628: if (status == GSL_SUCCESS) printf("Minimum found\n");
! 629: t_x11 = gsl_vector_get(s->x, 0);
! 630: t_x12 = gsl_vector_get(s->x, 1);
! 631: t_x22 = gsl_vector_get(s->x, 2);
! 632: t_y1 = gsl_vector_get(s->x, 3);
! 633: t_y2 = gsl_vector_get(s->x, 4);
! 634: printf("%d, %g, %g, %g, %g, %g, %g\n", iter, t_x11, t_x12, t_x22, t_y1, t_y2, s->f);
! 635: move_t_status = move_t(g_x11, g_x12, g_x22, g_y1, g_y2, t_x11, t_x12, t_x22, t_y1, t_y2, values);
! 636: if (move_t_status == MOVE_T_SUCCESS) {
! 637: g_x11 = t_x11;
! 638: g_x12 = t_x12;
! 639: g_x22 = t_x22;
! 640: g_y1 = t_y1;
! 641: g_y2 = t_y2;
! 642: } else {
! 643: printf("move_t failed\n");
! 644: g_x11 = move_t_points[0];
! 645: g_x12 = move_t_points[1];
! 646: g_x22 = move_t_points[2];
! 647: g_y1 = move_t_points[3];
! 648: g_y2 = move_t_points[4];
! 649: }
! 650: printf("points = [%g, %g, %g, %g, %g]\n", g_x11, g_x12, g_x22, g_y1, g_y2);
! 651: printf("values = [%g, %g, %g, %g]\n", values[0], values[1], values[2], values[3]);
! 652: printf("grad : "); gsl_vector_show(grad_v);
! 653: printf("norm(grad) : %f\n", gsl_blas_dnrm2(grad_v));
! 654: } while (status == GSL_CONTINUE && iter < SEARCH_MIN_ITER);
! 655: gsl_multimin_fdfminimizer_free(s);
! 656: gsl_vector_free(x);
! 657: }
! 658:
! 659: double my_f(const gsl_vector *v, void *params)
! 660: {
! 661: extern double values[N_VALUES];
! 662: extern double g_x11, g_x12, g_x22, g_y1, g_y2;
! 663: double o_x11, o_x12, o_x22, o_y1, o_y2;
! 664: double n_x11, n_x12, n_x22, n_y1, n_y2;
! 665: double val[N_VALUES];
! 666: int i;
! 667: for (i = 0; i < N_VALUES; i++)
! 668: val[i] = values[i];
! 669: o_x11 = g_x11;
! 670: o_x12 = g_x12;
! 671: o_x22 = g_x22;
! 672: o_y1 = g_y1;
! 673: o_y2 = g_y2;
! 674: n_x11 = gsl_vector_get(v, 0);
! 675: n_x12 = gsl_vector_get(v, 1);
! 676: n_x22 = gsl_vector_get(v, 2);
! 677: n_y1 = gsl_vector_get(v, 3);
! 678: n_y2 = gsl_vector_get(v, 4);
! 679: #ifdef DEBUG_RK
! 680: printf("call my_f([%g, %g, %g, %g, %g])\n", n_x11, n_x12, n_x22, n_y1, n_y2);
! 681: #endif
! 682: move_t(g_x11, g_x12, g_x22, g_y1, g_y2, n_x11, n_x12, n_x22, n_y1, n_y2, val);
! 683: g_x11 = o_x11;
! 684: g_x12 = o_x12;
! 685: g_x22 = o_x22;
! 686: g_y1 = o_y1;
! 687: g_y2 = o_y2;
! 688: return val[0];
! 689: }
! 690:
! 691: void my_df(const gsl_vector *v, void *param, gsl_vector *df)
! 692: {
! 693: extern double g_x11, g_x12, g_x22, g_y1, g_y2;
! 694: extern double values[N_VALUES];
! 695: extern gsl_vector *grad_v;
! 696:
! 697: double x11, x12, x22, y1, y2;
! 698: x11 = gsl_vector_get(v, 0);
! 699: x12 = gsl_vector_get(v, 1);
! 700: x22 = gsl_vector_get(v, 2);
! 701: y1 = gsl_vector_get(v, 3);
! 702: y2 = gsl_vector_get(v, 4);
! 703: #ifdef DEBUG_RK
! 704: printf("call my_df([%g, %g, %g, %g, %g])\n", x11, x12, x22, y1, y2);
! 705: #endif
! 706:
! 707: /* 点 v での関数値を計算 */
! 708: double o_x11, o_x12, o_x22, o_y1, o_y2;
! 709: double val[N_VALUES];
! 710: int i;
! 711:
! 712: /* move_all で更新されてしまうので、 g_* をコピー */
! 713: for (i = 0; i < N_VALUES; i++)
! 714: val[i] = values[i];
! 715: o_x11 =g_x11;
! 716: o_x12 =g_x12;
! 717: o_x22 =g_x22;
! 718: o_y1 = g_y1;
! 719: o_y2 = g_y2;
! 720:
! 721: move_t(g_x11, g_x12, g_x22, g_y1, g_y2, gsl_vector_get(v, 0), gsl_vector_get(v, 1), gsl_vector_get(v, 2), gsl_vector_get(v, 3), gsl_vector_get(v, 4), val);
! 722: /* g_* を元の値に戻す */
! 723: g_x11 =o_x11;
! 724: g_x12 =o_x12;
! 725: g_x22 =o_x22;
! 726: g_y1 = o_y1;
! 727: g_y2 = o_y2;
! 728: grad(x11, x12, x22, y1, y2, 1, val);
! 729: for (i = 0; i < DIM; i++)
! 730: gsl_vector_set(df, i, gsl_vector_get(grad_v, i));
! 731: }
! 732:
! 733: void my_fdf(const gsl_vector *x, void *params, double *f, gsl_vector *df)
! 734: {
! 735: *f = my_f(x, params);
! 736: my_df(x, params, df);
! 737: }
! 738:
! 739: void grad(double x11, double x12, double x22, double y1, double y2, double r, double *val)
! 740: {
! 741: extern gsl_vector *grad_v;
! 742: extern gsl_matrix *pf11_m, *pf12_m, *pf22_m, *pf1_m, *pf2_m;
! 743: int i;
! 744: double t;
! 745:
! 746: for (i = 0, t = 0; i < N_VALUES; i++)
! 747: t += gsl_matrix_get(pf11_m, 0, i) * val[i];
! 748: gsl_vector_set(grad_v, 0, t);
! 749:
! 750: for (i = 0, t = 0; i < N_VALUES; i++)
! 751: t += gsl_matrix_get(pf12_m, 0, i) * val[i];
! 752: gsl_vector_set(grad_v, 1, t);
! 753:
! 754: for (i = 0, t = 0; i < N_VALUES; i++)
! 755: t += gsl_matrix_get(pf22_m, 0, i) * val[i];
! 756: gsl_vector_set(grad_v, 2, t);
! 757:
! 758: for (i = 0, t = 0; i < N_VALUES; i++)
! 759: t += gsl_matrix_get(pf1_m, 0, i) * val[i];
! 760: gsl_vector_set(grad_v, 3, t);
! 761:
! 762: for (i = 0, t = 0; i < N_VALUES; i++)
! 763: t += gsl_matrix_get(pf2_m, 0, i) * val[i];
! 764: gsl_vector_set(grad_v, 4, t);
! 765:
! 766: }
! 767:
! 768: void gsl_vector_show(gsl_vector *v)
! 769: {
! 770: int i;
! 771: for (i = 0; i < v->size; i++)
! 772: printf("%f ", gsl_vector_get(v, i));
! 773: printf("\n");
! 774: }
! 775:
! 776: void show_v(double *v, int n)
! 777: {
! 778: int i;
! 779: printf("[");
! 780: for (i = 0; i < n; i++) {
! 781: printf("%g", v[i]);
! 782: if (i < n - 1)
! 783: printf(", ");
! 784: }
! 785: printf("]\n");
! 786: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to test1.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / packages / doc / nk_fb_gen_c