とりあえず,MathJaxのテスト。デフォルトをDisplay modeにして
\[ y(s) = f(s) + \int^b_a g(s, t) y(t) dt \]
数式が出るかしら? ついでにインライン表示も・・・まぁこんなもんか。
で,コード例も一つ。$n$次正方行列と次元ベクトルから作られる連立一次方程式
\[ A\mathbf{x} = \mathbf{b} \]
を未知ベクトルについて解くLAPACKプログラム。DGESVを使う最も簡単なもの。解の表示をコメントにしてあるのはでかい次元数のベンチマーク用に使ったからです。適宜はずしてお使い下さいませ。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 |
/*************************************************/ /* LAPACK Tutorial */ /* Sample Program of DGESV */ /* with Intel Math Kernel */ /* Last Update: 2013-01-21 (Mon) T.Kouya */ // icc -I/opt/intel/include -I/opt/intel/mkl/include -L/opt/intel/mkl/lib/intel64 -L/opt/intel/lib/intel64 blas1_mkl.c -Wl,--start-group /opt/intel/mkl/lib/intel64/libmkl_intel_lp64.a /opt/intel/mkl/lib/intel64/libmkl_sequential.a /opt/intel/mkl/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread /*************************************************/ #include <stdio.h> #include <stdlib.h> #include "mkl.h" // for Intel Math Kernel Library #include "mkl_cblas.h" #include "mkl_lapacke.h" #define min(i, j) (((i) < (j)) ? (i) : (j)) #define max(i, j) (((i) > (j)) ? (i) : (j)) int main() { MKL_INT i, j, dim; // dimension of vectors MKL_INT *pivot, info; double *ma, *vx, *vb; double alpha, beta; char str_mkl_version[1024]; double running_time; // print MKL version MKL_Get_Version_String(str_mkl_version, 1024); printf("%s\n", str_mkl_version); // input dimension printf("Dim = "); scanf("%d", &dim); if(dim <= 0) { printf("Illigal dimension! (dim = %d)\n", dim); return EXIT_FAILURE; } // Initialize pivot = (MKL_INT *)calloc(sizeof(MKL_INT), dim); ma = (double *)calloc(sizeof(double), dim * dim); vx = (double *)calloc(sizeof(double), dim); vb = (double *)calloc(sizeof(double), dim); // input va and vb for(i = 0; i < dim; i++) { for(j = 0; j < dim; j++) { // A = hirbert matrix //*(ma + i * dim + j) = 1.0 / (double)(i + j + 1); // A = frank matrix *(ma + i * dim + j) = (double)(dim - max(i, j)); } // x = answer *(vx + i) = (double) (i + 1); } // b := A * x alpha = 1.0; beta = 0.0; cblas_dgemv(CblasRowMajor, CblasNoTrans, dim, dim, alpha, ma, dim, vx, 1, beta, vb, 1); // print /* printf("A = \n"); for(i = 0; i < dim; i++) { printf("%3d: ", i); for(j = 0; j < dim; j++) printf("%10f ", *(ma + i * dim + j)); printf("\n"); } // print printf("x = \n"); for(i = 0; i < dim; i++) { printf("%3d: ", i); printf("%10f ", *(vx + i)); printf("\n"); } // print printf("b = \n"); for(i = 0; i < dim; i++) { printf("%3d: ", i); printf("%10f ", *(vb + i)); printf("\n"); } */ // solve A * X = C -> C := X running_time = dsecnd(); // start info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, dim, 1, ma, dim, pivot, vb, 1); running_time = dsecnd() - running_time; // end // error occurs if info > 0 if(info > 0) { printf("pivot ma[%d][%d] is zero!\n", info, info); return EXIT_FAILURE; } printf("Dimension: %d, Running Time(s): %g\n", dim, running_time); // print /* printf("calculated x = \n"); for(i = 0; i < dim; i++) { printf("%3d: ", i); printf("%10f ", *(vb + i)); printf("\n"); } // diff printf("x - calculated x = \n"); for(i = 0; i < dim; i++) { printf("%3d: ", i); printf("%10f ", fabs(*(vx + i) - *(vb + i))); printf("\n"); } */ // free free(pivot); free(ma); free(vb); free(vx); return EXIT_SUCCESS; } |