• xunhuan Jacobi method SVD

    单纯的cpp程序 ，svd， 循环Jacobimethod

    1. #include <stdio.h>
    2. #include <math.h>
    3. #include <stdlib.h>
    4.  
    5. void printFloatMatrixRowMajor(float* A, int M, int N, int lda);
    6. void printFloatMatrixColMajor(float* A, int M, int N, int lda);
    7.  
    8.  
    9.  
    10. int jcbi_roll(float* a, int n, float* v, float eps, int jt) {
    11.         int i, j, p, q, u, w, t, s, loops;
    12.         float fm, cn, sn, omega, x, y, d;
    13.         loops = 1;
    14.         int contiCount = 0;
    15.  
    16.         for (i = 0; i <= n - 1; i++)//Let V = I;
    17.         {
    18.             v[i * n + i] = 1.0;
    19.             for (j = 0; j <= n - 1; j++)
    20.                 if (i != j) v[i * n + j] = 0.0;
    21.         }
    22.  
    23.         while (1 == 1)
    24.         {
    25.             fm = 0.0;
    26.             for (p = 1; p <= n - 1; p++)
    27.             {
    28.                 for (q = 0; q <= p - 1; q++)
    29.                 {
    30.                     if (loops > jt) { printf("\n\ncontiCount = %d\n", contiCount);  return(-1); }// 超过迭代次数，依然没有收敛，退出，返回-1
    31.  
    32.                     loops = loops + 1;//loopCount++ ;
    33.                     u = p * n + q; //A(p, q)
    34.                     w = p * n + p; //A(p, p)
    35.                     t = q * n + p; //A(q, p)
    36.                     s = q * n + q; //A(q, q)
    37.                     if (fabs(a[u]) < eps) { contiCount++; continue; }
    38.  
    39.                     //计算sin(2θ)
    40.                     x = -a[u];                          // m = -A(p,q)
    41.                     y = (a[s] - a[w]) / 2.0;            // n = 1/2(A(q,q)-A(pp))
    42.                     omega = x / sqrt(x * x + y * y);    // m/(sqrt(m^2 + n^2))             
    43.                     if (y < 0.0) omega = -omega;        //sgn(n)*              sin(2θ)=omega
    44.                     //计算sin(θ)
    45.                     sn = 1.0 + sqrt(1.0 - omega * omega);
    46.                     sn = omega / sqrt(2.0 * sn);        // sin(θ)
    47.                     //计算cos(θ)
    48.                     cn = sqrt(1.0 - sn * sn);           // cos(θ)
    49.                     fm = a[w];                          // A(p,p)
    50.                     //计算 四个元素
    51.                     a[w] = fm * cn * cn + a[s] * sn * sn + a[u] * omega; // A(p,p)' = A(p,p)cos(θ)^2 + A(q,q)sin(θ)^2 + A(p,q)sin(2θ)
    52.                     a[s] = fm * sn * sn + a[s] * cn * cn - a[u] * omega; // A(q,q)' = A(p,p)sin(θ)^2 + A(q,q)cos(θ)^2 - A(p,q)sin(2θ)
    53.                     a[u] = 0.0;                                          // A(p,q)' = 0;
    54.                     a[t] = 0.0;                                          // A(q,p)' = 0;
    55.                     //计算第p行 和第q行
    56.                     for (j = 0; j <= n - 1; j++)            // updat row p  and  row q;
    57.                         if ((j != p) && (j != q))
    58.                         {
    59.                             u = p * n + j;               // 遍历第 p 行的每一个元素
    60.                             w = q * n + j;               // 遍历第 q 行的每一个元素
    61.                             fm = a[u];                    // 
    62.                             a[u] = fm * cn + a[w] * sn;    // 自己的上次值 累加上 伴随行的元素上次值
    63.                             a[w] = -fm * sn + a[w] * cn;    // 自己的原值 累加上 伴随行的元素值
    64.                         }
    65.                     //计算第p列和第q列
    66.                     for (i = 0; i <= n - 1; i++)            // update column p  and  column q
    67.                         if ((i != p) && (i != q))
    68.                         {
    69.                             u = i * n + p;               // A(i, p) 遍历 p-th column 
    70.                             w = i * n + q;               // A(i, q) 遍历 q-th column
    71.                             fm = a[u];                    //
    72.                             a[u] = fm * cn + a[w] * sn;     // 自己的上次值 累加上 伴随行的元素上次值
    73.                             a[w] = -fm * sn + a[w] * cn;    // 自己的上次值 累加上 伴随行的元素上次值
    74.                         }
    75.                     for (i = 0; i <= n - 1; i++)            // update EigenVectors p-th  and q-th
    76.                     {
    77.                         u = i * n + p;                      // V(i, p) 遍历 p-th column
    78.                         w = i * n + q;                      // V(i, q) 遍历 q-th column
    79.                         fm = v[u];                          //
    80.                         v[u] = fm * cn + v[w] * sn;         // 
    81.                         v[w] = -fm * sn + v[w] * cn;        // 
    82.                     } 
    83.                 }///LLLL
    84.             }
    85.         }
    86.  
    87.         printf("\n\ncontiCount=%d\n", contiCount);
    88.  
    89.         return(1);
    90.     }
    91.  
    92.  
    93. int jcbi(float* a, int n, float* v, float eps, int jt){
    94.     int i, j, p, q, u, w, t, s, loops;
    95.     float fm, cn, sn, omega, x, y, d;
    96.     loops = 1;
    97.  
    98.     for (i = 0; i <= n - 1; i++)//Let V = I;
    99.     {
    100.         v[i * n + i] = 1.0;
    101.         for (j = 0; j <= n - 1; j++)
    102.             if (i != j) v[i * n + j] = 0.0;
    103.     }
    104.  
    105.     while (1 == 1)
    106.     {
    107.         fm = 0.0;
    108.         for (i = 1; i <= n - 1; i++) {//only check the lower triangle matrix of D_i, without the diag elements
    109.             for (j = 0; j <= i - 1; j++)// q 恒小于 p;   q 不可能等于 p
    110.             {
    111.                 d = fabs(a[i * n + j]);
    112.                 if ((i != j) && (d > fm))
    113.                 {
    114.                     fm = d;//当前最大值
    115.                     p  = i;//行号
    116.                     q  = j;//列号
    117.                 }
    118.             }
    119.         }
    120.  
    121.         if (fm < eps)  return(1);//the max element is less than eps, tolerance
    122.         if ( loops > jt)  return(-1);// 超过迭代次数，依然没有收敛，退出，返回-1
    123.  
    124.         loops = loops + 1;//loopCount++ ;
    125.         u = p * n + q; //A(p, q)
    126.         w = p * n + p; //A(p, p)
    127.         t = q * n + p; //A(q, p)
    128.         s = q * n + q; //A(q, q)
    129.  
    130. //计算sin(2θ)
    131.         x = -a[u];                          // m = -A(p,q)
    132.         y = (a[s] - a[w]) / 2.0;            // n = 1/2(A(q,q)-A(pp))
    133.         omega = x / sqrt(x * x + y * y);    // m/(sqrt(m^2 + n^2))             
    134.         if (y < 0.0) omega = -omega;        //sgn(n)*              sin(2θ)=omega
    135. //计算sin(θ)
    136.         sn = 1.0 + sqrt(1.0 - omega * omega);
    137.         sn = omega / sqrt(2.0 * sn);        // sin(θ)
    138. //计算cos(θ)
    139.         cn = sqrt(1.0 - sn * sn);           // cos(θ)
    140.         fm = a[w];                          // A(p,p)
    141. //计算 四个元素
    142.         a[w] = fm * cn * cn + a[s] * sn * sn + a[u] * omega; // A(p,p)' = A(p,p)cos(θ)^2 + A(q,q)sin(θ)^2 + A(p,q)sin(2θ)
    143.         a[s] = fm * sn * sn + a[s] * cn * cn - a[u] * omega; // A(q,q)' = A(p,p)sin(θ)^2 + A(q,q)cos(θ)^2 - A(p,q)sin(2θ)
    144.         a[u] = 0.0;                                          // A(p,q)' = 0;
    145.         a[t] = 0.0;                                          // A(q,p)' = 0;
    146. //计算第p行 和第q行
    147.         for (j = 0; j <= n - 1; j++)            // updat row p  and  row q;
    148.             if ((j != p) && (j != q))
    149.             {
    150.                 u    = p * n + j;               // 遍历第 p 行的每一个元素
    151.                 w    = q * n + j;               // 遍历第 q 行的每一个元素
    152.                 fm   = a[u];                    // 
    153.                 a[u] =  fm * cn + a[w] * sn;    // 自己的上次值 累加上 伴随行的元素上次值
    154.                 a[w] = -fm * sn + a[w] * cn;    // 自己的原值 累加上 伴随行的元素值
    155.             }
    156. //计算第p列和第q列
    157.         for (i = 0; i <= n - 1; i++)            // update column p  and  column q
    158.             if ((i != p) && (i != q))
    159.             {
    160.                 u    = i * n + p;               // A(i, p) 遍历 p-th column 
    161.                 w    = i * n + q;               // A(i, q) 遍历 q-th column
    162.                 fm   = a[u];                    //
    163.                 a[u] = fm * cn + a[w] * sn;     // 自己的上次值 累加上 伴随行的元素上次值
    164.                 a[w] = -fm * sn + a[w] * cn;    // 自己的上次值 累加上 伴随行的元素上次值
    165.             }
    166.         for (i = 0; i <= n - 1; i++)            // update EigenVectors p-th  and q-th
    167.         {
    168.             u = i * n + p;                      // V(i, p) 遍历 p-th column
    169.             w = i * n + q;                      // V(i, q) 遍历 q-th column
    170.             fm = v[u];                          //
    171.             v[u] = fm * cn + v[w] * sn;         // 
    172.             v[w] = -fm * sn + v[w] * cn;        // 
    173.         }
    174.     }
    175.  
    176.     return(1);
    177. }
    178.  
    179. void transpositionColMajor(float* A, int M, int N, float* At, int lda, int ldat) {
    180.     for (int i = 0; i < N; i++) {
    181.         for (int j = 0; j < M; j++) {
    182.             At[i + j * ldat] = A[j + i * lda];
    183.         }
    184.     }
    185. }
    186. void transpositionRowMajor(float* A, int M, int N, float* At, int lda, int ldat) {
    187.     for (int i = 0; i < N; i++) {
    188.         for (int j = 0; j < M; j++) {
    189.             At[i * ldat + j] = A[j * lda + i];
    190.         }
    191.     }
    192. }
    193. void SGemmColMajor(float* A, float* B, int M, int N, int K, float* C, int lda, int ldb, int ldc) {
    194.     for (int i = 0; i < M; i++) {
    195.         for (int j = 0; j < N; j++) {
    196.             C[i + j * ldc] = 0.0f;
    197.             for (int k = 0; k < K; k++) {
    198.                 C[i + j * ldc] += A[i + k * lda] * B[k + j * ldb];
    199.             }
    200.         }
    201.     }
    202. }
    203. void SGemmRowMajor(float* A, float* B, int M, int N, int K, float* C, int lda, int ldb, int ldc) {
    204.     for (int i = 0; i < M; i++) {
    205.         for (int j = 0; j < N; j++) {
    206.             C[i * ldc + j] = 0.0f;
    207.             for (int k = 0; k < K; k++) {
    208.                 C[i * ldc + j] += A[i * lda + k] * B[k * ldb + j];
    209.             }
    210.         }
    211.     }
    212. }
    213. void printFloatMatrixColMajor(float* A, int M, int N, int lda) {
    214.     for (int i = 0; i < M; i++) {
    215.         for (int j = 0; j < N; j++) {
    216.             printf(" %9.7f", A[i + j * lda]);
    217.         }
    218.         printf("\n\n");
    219.     }
    220. }
    221. void printFloatMatrixRowMajor(float* A, int M, int N, int lda) {
    222.     for (int i = 0; i < M; i++) {
    223.         for (int j = 0; j < N; j++) {
    224.             printf(" %6.4f", A[i * lda + j]);
    225.         }
    226.         printf("\n\n");
    227.     }
    228. }
    229. void gemmTest() {
    230.     int m, n, k;
    231.     m = 4;
    232.     n = 2;
    233.     k = 3;
    234.     float A[12] = { 5,3,6,0,2,8,0,1,4,2,4,6 };
    235.     float B[6] = { 2,1,3,4,3,2 };
    236.     float C[8] = { 0.0 };
    237.     printf("\nA=\n");
    238.     printFloatMatrixColMajor(A, m, k, m);
    239.     printf("\nB=\n");
    240.     printFloatMatrixColMajor(B, k, n, k);
    241.     SGemmColMajor(A, B, m, n, k, C, m, k, m);
    242.     printf("\nC=\n");
    243.     printFloatMatrixColMajor(C, m, n, m);
    244. }
    245. void transpositionTest() {
    246.     int m, n, k;
    247.     m = 4;
    248.     n = 3;
    249.  
    250.     float  A[12] = { 5,3,6,0,2,8,0,1,4,2,4,6 };
    251.     float At[12];
    252.     //void transpositionColMajor(float* A, int M, int N, float* At, int lda, int ldat)
    253.     transpositionColMajor(A, m, n, At, m, n);
    254.  
    255.     printf("\nA=\n");
    256.     printFloatMatrixColMajor(A, m, n, m);
    257.     printf("\nAt=\n");
    258.     printFloatMatrixColMajor(At, n, m, n);
    259. }
    260. void svd_sym_Little_Test() {
    261.     int n = 6;
    262.  
    263.     //float  A[36] = { 2, -1, 0, -1, 2, -1, 0, -1, 2 };
    264.     float A[36] = {
    265.    1.9757,   2.4581,   1.4726, 1.3031, 1.5231, 0.8238,
    266.    2.4581,   3.5969,   1.9414, 2.3095, 2.3901, 1.4690,
    267.    1.4726,   1.9414,   1.5127, 1.2891, 0.9270, 0.5674,
    268.    1.3031,   2.3095,   1.2891, 1.8225, 1.4753, 1.0712,
    269.    1.5231,   2.3901,   0.9270, 1.4753, 2.0327, 0.9849,
    270.    0.8238,   1.4690,   0.5674, 1.0712, 0.9849, 0.9280
    271.     };
    272.  
    273.  
    274.     printf("\nA=\n");
    275.     printFloatMatrixRowMajor(A, n, n, n);
    276.  
    277.     int   jt;
    278.     float eps;
    279.     int result = 0;
    280.     float eigenVs[36];
    281.  
    282.     jt = 100;
    283.     eps = 0.000001;
    284.     result = jcbi(A, n, eigenVs, eps, jt);
    285.  
    286.     printf("\neigenValus=\n");
    287.     printFloatMatrixRowMajor(A, n, n, n);
    288.  
    289.  
    290.     printf("\neigenVs=\n");
    291.     printFloatMatrixRowMajor(eigenVs, n, n, n);
    292. }
    293. void svd_Little_mGEn_Test() {
    294.     int m, n, k;
    295.     m = 4;
    296.     n = 3;
    297.  
    298.     //float  A[12] = { 5,3,6,0,2,8,0,1,4,2,4,6 };
    299.     float  A[12] = { 5,2,4,3,8,2,6,0,4,0,1,6 };
    300.  
    301.     float At[12];
    302.     //void transpositionRowMajor(float* A, int M, int N, float* At, int lda, int ldat)
    303.     transpositionRowMajor(A, m, n, At, n, m);
    304.  
    305.     printf("\nA=\n");
    306.     printFloatMatrixRowMajor(A, m, n, n);
    307.     printf("\nAt=\n");
    308.     printFloatMatrixRowMajor(At, n, m, m);
    309.  
    310.     float A_x_At[16];
    311.     float At_x_A[9];
    312.  
    313.     //void SGemmColMajor(float* A, float* B, int M, int N, int K, float* C, int lda, int ldb, int ldc)
    314.     SGemmRowMajor(At, A, n, n, m, At_x_A, m, n, n);
    315.     printf("\nAt_x_A=\n");
    316.     printFloatMatrixRowMajor(At_x_A, n, n, n);
    317.  
    318.     int jt;
    319.     float eps;
    320.     int result = 0;
    321.     float v_At_x_A[9];
    322.  
    323.     jt = 100;
    324.     eps = 0.000001;
    325.  
    326.     result = jcbi(At_x_A, n, v_At_x_A, eps, jt);
    327.     printf("\nSigma in At_x_A=\n");
    328.     printFloatMatrixRowMajor(At_x_A, n, n, n);
    329.     printf("\nV of At_x_A=\n");
    330.     printFloatMatrixRowMajor(v_At_x_A, n, n, n);
    331.  
    332. }
    333. void svd_Little_mLTn_Test() {
    334.     int m, n, k;
    335.     m = 4;
    336.     n = 3;
    337.  
    338.     //float  A[12] = { 5,3,6,0,2,8,0,1,4,2,4,6 };
    339.     float  A[12] = { 5,2,4,3,8,2,6,0,4,0,1,6 };
    340.  
    341.     float At[12];
    342.     //void transpositionColMajor(float* A, int M, int N, float* At, int lda, int ldat)
    343.     transpositionColMajor(A, m, n, At, m, n);
    344.  
    345.     printf("\nA=\n");
    346.     printFloatMatrixColMajor(A, m, n, m);
    347.     printf("\nAt=\n");
    348.     printFloatMatrixColMajor(At, n, m, n);
    349.  
    350.     float A_x_At[16];
    351.     float At_x_A[9];
    352.  
    353.     //void SGemmColMajor(float* A, float* B, int M, int N, int K, float* C, int lda, int ldb, int ldc)
    354.     /****
    355.     SGemmColMajor(A, At, m, m, n, A_x_At, m, n, m);
    356.     printf("\nA_x_At=\n");
    357.     printFloatMatrixColMajor(A_x_At, m, m, m);
    358.     *****/
    359.     SGemmColMajor(At, A, n, n, m, At_x_A, n, m, n);
    360.     printf("\nAt_x_A=\n");
    361.     printFloatMatrixColMajor(At_x_A, n, n, n);
    362.  
    363.     int jt;
    364.     float eps;
    365.     int result = 0;
    366.     float v_A_x_At[16];
    367.     float v_At_x_A[9];
    368.  
    369.     jt = 100;
    370.     eps = 0.000001;
    371.  
    372.     /****
    373.     result = jcbi(A_x_At, m, v_A_x_At, eps, jt);
    374.  
    375.     printf("\nA_x_At=\n");
    376.     printFloatMatrixColMajor(A_x_At, m, m, m);
    377.     printf("\nV of A_x_At=\n");
    378.     printFloatMatrixColMajor(v_A_x_At, m, m, m);
    379.     *****/
    381.     result = jcbi(At_x_A, n, v_At_x_A, eps, jt);
    382.     printf("\nAt_x_A=\n");
    383.     printFloatMatrixColMajor(At_x_A, n, n, n);
    384.     printf("\nV of At_x_A=\n");
    385.     printFloatMatrixColMajor(v_At_x_A, n, n, n);
    386.  
    387. }
    388. void printFloatVector(float* X, int n) {
    389.     for (int idx = 0; idx < n; idx++) {
    390.         printf(" %9.7e,", X[idx]);
    391.     }
    392.     printf("\n");
    393. }
    394. void svd_Big_Test() {
    395.     int m, n, k;
    396.     int batchSize = 1;
    397.     m = 300;
    398.     n = 181;
    399.  
    400.     float* A;
    401.     float* At;
    402.     A = (float*)malloc(sizeof(float) * m * n * batchSize);
    403.     At = (float*)malloc(sizeof(float) * n * m * batchSize);
    404.     const char* fileName = "test_svd.bin";
    405.     //void transpositionColMajor(float* A, int M, int N, float* At, int lda, int ldat)
    406.     //transpositionColMajor(A, m, n, At, m, n);
    407.     FILE* fp = fopen(fileName, "rb");
    408.  
    409.     if (fp == NULL) {
    410.         printf("FILE OPEN FAIL!\n");
    411.         return;
    412.     }
    413.  
    414.     fread(A, sizeof(float), m * n * batchSize, fp);
    415.     //fread(A, sizeof(float), 7, fp);
    416.  
    417.     transpositionColMajor(A, m, n, At, m, n);
    418.  
    419.     printf("A=\n\n");
    420.     printFloatVector(A, 7);
    421.  
    422.  
    423.     printf("\nA=\n");
    424.     //printFloatMatrixColMajor(A, m, n, m);
    425.     printf("\nAt=\n");
    426.     //printFloatMatrixColMajor(At, n, m, n);
    427.  
    428.     float* A_x_At;
    429.     float* At_x_A;
    430.  
    431.     A_x_At = (float*)malloc(m * m * sizeof(float) * batchSize);
    432.     At_x_A = (float*)malloc(n * n * sizeof(float) * batchSize);
    433.  
    434.     //void SGemmColMajor(float* A, float* B, int M, int N, int K, float* C, int lda, int ldb, int ldc)
    435.     SGemmColMajor(A, At, m, m, n, A_x_At, m, n, m);
    436.     printf("\nA_x_At=\n");
    437.     printFloatMatrixColMajor(A_x_At, m, m, m);
    438.  
    439.     SGemmColMajor(At, A, n, n, m, At_x_A, n, m, n);
    440.     printf("\nAt_x_A=\n");
    441.     printFloatMatrixColMajor(At_x_A, n, n, n);
    442.  
    443.     int jt;
    444.     float eps;
    445.     int result = 0;
    446.     float* v_A_x_At;
    447.     float* v_At_x_A;
    448.     v_A_x_At = (float*)malloc(m * m * sizeof(float) * batchSize);
    449.     v_At_x_A = (float*)malloc(n * n * sizeof(float) * batchSize);
    450.  
    451.     jt = 100;
    452.     eps = 0.000001;
    453.     result = jcbi(A_x_At, m, v_A_x_At, eps, jt);
    454.  
    455.     printf("\nA_x_At=\n");
    456.     printFloatMatrixColMajor(A_x_At, m, m, m);
    457.     printf("\nU of A_x_At=\n");
    458.     printFloatMatrixColMajor(v_A_x_At, m, m, m);
    459.  
    460.  
    461.     result = jcbi(At_x_A, n, v_At_x_A, eps, jt);
    462.     printf("\nAt_x_A=\n");
    463.     printFloatMatrixColMajor(At_x_A, n, n, n);
    464.     printf("\nV of At_x_A=\n");
    465.     printFloatMatrixColMajor(v_At_x_A, n, n, n);
    466.     /****/
    467. }
    468. void svd_sym_Test() {
    469.     int m, n;
    470.     m = 7;
    471.     n = 7;
    472.     float A[49] = {
    473.                                   6, 6, 0, 9, 7, 2, 6,
    474.                                   6, 5, 9, 4, 2, 5, 5,
    475.                                   0, 9, 8, 3, 5, 0, 6,
    476.                                   9, 4, 3, 8, 9, 9, 4,
    477.                                   7, 2, 5, 9, 0, 6, 0,
    478.                                   2, 5, 0, 9, 6, 4, 4,
    479.                                   6, 5, 6, 4, 0, 4, 9
    480.     };
    481.  
    482.  
    483.     /**
    484.     float A[4096] = {
    485.    9, 10,  7,  9,  2,  0, 19, 12, 16,  4,  0,  1, 12,  4,  7,  4, 14, 18, 11,  1, 12, 19, 14, 12, 12,  0,  9, 10,  3,  8,  6, 16, 15,  5, 10,  4, 10, 14,  4, 10, 13, 10,  4,  1,  7,  3,  3, 11,  1, 15,  7, 19, 19, 14,  0, 13,  4, 19, 19,  7,  6,  3,  0, 12
    486. , 10, 14, 11, 10,  0, 11, 19, 15, 13, 19,  5, 14,  1,  3,  8,  7,  9, 14, 12,  2, 16, 14,  3, 19,  1,  2,  3, 18, 19,  9, 15,  0,  2, 18, 18, 15,  2, 14, 11, 12, 16, 11, 10,  6, 11,  5,  9,  5,  9,  1, 14, 17,  4,  4, 10,  0, 17,  7,  1, 18, 16,  5,  3, 14
    487. ,  7, 11, 16, 19, 13,  1,  2,  1, 13,  4, 15, 17,  4,  0,  1, 18, 13,  4, 14, 13, 10, 12, 12, 19,  8, 10,  0,  4,  4,  1, 17, 14,  2, 14,  9, 14, 15,  2, 18, 13, 16,  2, 15,  7, 14,  4,  1,  0,  2,  8,  1,  1,  3,  6, 17,  9, 15, 12,  4,  4, 13,  8, 14,  1
    488. ,  9, 10, 19,  0,  2,  4, 11, 16,  6,  6, 18,  6, 16,  2,  9, 11,  5,  4,  5,  1, 11,  9, 13,  3,  2,  3,  3, 16,  6,  9, 11,  0, 19, 11, 10, 13,  6,  6, 13,  3, 12,  6, 18,  4,  8,  9, 13, 14, 13, 16, 19, 14, 12, 11, 18, 11,  9,  5,  4,  9,  3,  4,  4,  1
    489. ,  2,  0, 13,  2, 10, 11, 16, 19, 17, 17, 11, 15, 15, 19, 19, 10,  2, 13,  1,  9,  2,  6,  4,  9, 17,  9, 17,  5, 12, 16,  9,  5, 18,  7,  6, 12,  9, 18, 13,  8, 12,  9, 10,  0,  8, 12, 14, 17,  0,  5, 14,  8, 11,  9,  0,  4,  2,  5,  0,  5,  1,  0,  7,  8
    490. ,  0, 11,  1,  4, 11, 10, 11, 15,  8, 15, 17,  6,  9,  3, 16,  4, 16,  0, 15, 16,  4,  0, 18,  3,  9, 14,  6,  7, 11,  1, 11, 17,  4, 15, 10, 13,  6, 15, 14, 12,  5,  3,  6,  1,  4,  7,  9, 16, 12,  2, 17,  8,  7,  5, 17,  6,  3, 14,  1, 16,  4, 13, 10, 10
    491. , 19, 19,  2, 11, 16, 11,  5,  5,  7,  6, 17, 12,  5, 14, 16, 18, 19,  6, 16, 16, 19,  9,  1,  0,  9,  6, 15,  7, 18,  6, 18,  4,  8, 12, 16, 15, 13,  7, 11,  8, 12,  3,  1, 19,  8, 13,  8,  6,  6,  9, 11, 11, 11,  7, 16,  2, 19, 14, 16, 19,  5,  7,  6,  8
    492. , 12, 15,  1, 16, 19, 15,  5,  2,  3, 17, 10,  0,  9,  7,  6,  9,  1,  6,  7, 16, 16, 17, 17,  9, 18,  4,  2,  9, 16,  1, 16,  4,  0,  5, 18,  6,  0,  6, 12, 14,  4,  8, 19,  3,  7,  5,  7, 11, 12, 15, 13,  8,  6, 15, 11, 18, 11,  7,  6,  3, 17, 16, 17,  4
    493. , 16, 13, 13,  6, 17,  8,  7,  3,  0,  7, 13,  2,  8, 12, 19,  5, 14,  1,  4, 17,  0,  5,  8, 13,  6,  0, 17,  7,  8,  5,  4, 14, 12,  2, 14,  1,  1,  0,  5, 13, 10,  9,  6,  1, 17,  8,  5,  6, 10, 14,  3,  5,  6,  0,  7,  4,  8, 17, 19, 10,  6,  2,  7, 12
    494. ,  4, 19,  4,  6, 17, 15,  6, 17,  7, 12,  4,  1, 12,  1,  9,  5, 10,  0, 16,  1,  3,  0,  0, 18,  7,  4,  7, 11,  0, 14,  2, 15,  2,  2, 11,  9, 11, 16,  8, 16,  9,  7, 17,  9,  4,  3,  1,  5,  7, 19,  5,  5,  9, 17, 14, 10, 15, 19,  5,  9, 15,  0, 10, 19
    495. ,  0,  5, 15, 18, 11, 17, 17, 10, 13,  4,  3, 10,  2, 17, 19,  6,  5,  5,  4, 18,  8, 15, 14, 18,  9,  2,  2, 18, 10, 17, 14,  7, 13,  1, 11,  6, 17,  6,  1, 18,  9,  4, 15, 11,  4, 11, 16, 18,  3,  9, 18,  2,  2,  1,  4, 17,  4, 16, 14,  5,  6, 19, 15, 18
    496. ,  1, 14, 17,  6, 15,  6, 12,  0,  2,  1, 10,  4, 19,  5,  6, 19,  4,  1,  0,  2,  1,  9,  2, 12,  5, 17,  6,  9, 14, 14,  2,  4, 17,  4,  2,  6, 10, 18,  5,  8, 10, 14, 18, 13,  8,  7,  1,  1,  6, 14, 17,  4, 19, 12,  4, 15, 15, 16,  3, 11,  9,  3, 10,  1
    497. , 12,  1,  4, 16, 15,  9,  5,  9,  8, 12,  2, 19,  9,  6, 11,  7, 15,  5,  7, 19,  5,  1, 12,  8, 13,  3, 14, 11, 16, 10, 10, 13, 17,  9, 11, 19, 12,  4,  8, 18, 10, 10,  0,  3,  5,  0,  0, 19,  5,  5, 16, 19, 18,  4,  2, 17,  8,  1,  3, 17, 10,  7,  0,  2
    498. ,  4,  3,  0,  2, 19,  3, 14,  7, 12,  1, 17,  5,  6,  4,  7,  2,  6,  2,  7,  1, 16, 10, 14,  1,  1,  6, 17,  7, 19,  8, 10,  0,  0, 10, 13,  5,  0,  0, 13, 18, 13, 14,  1, 11, 12,  8, 16,  8,  2, 17,  3, 13, 19,  1,  6,  0,  7, 12,  8,  6, 15,  7,  9, 19
    499. ,  7,  8,  1,  9, 19, 16, 16,  6, 19,  9, 19,  6, 11,  7, 16, 11, 14,  0,  7, 15,  7,  4,  1,  9,  5, 18,  8,  1,  0,  0, 17,  1,  1,  3, 12,  1, 17, 19, 14, 11, 10,  0, 17,  0, 11,  1, 14, 19, 17,  0,  7, 11, 16,  5, 10,  4,  2,  1, 10, 14, 15,  3,  7, 12
    500. ,  4,  7, 18, 11, 10,  4, 18,  9,  5,  5,  6, 19,  7,  2, 11, 14, 15,  7,  8, 10, 15, 10,  0, 10, 14, 16, 14,  4, 13, 13,  9, 19,  4, 15, 18, 12, 14,  5,  3,  0, 13, 12, 14,  3,  6, 19,  0,  6, 14, 13,  8, 12,  2,  7, 13, 16,  0,  0,  5, 16,  2,  8,  9, 11
    501. , 14,  9, 13,  5,  2, 16, 19,  1, 14, 10,  5,  4, 15,  6, 14, 15,  7,  3, 14,  2, 10, 15, 16,  0,  7,  8,  4, 17, 10,  3, 19, 15,  0,  7,  2, 14,  1,  4,  8, 18, 14,  4,  5, 10,  4, 17, 17,  6, 11,  4,  7,  3,  7,  9,  6, 15,  6, 10, 19, 19, 19,  7,  6,  5
    502. , 18, 14,  4,  4, 13,  0,  6,  6,  1,  0,  5,  1,  5,  2,  0,  7,  3,  2,  3,  0,  5,  2, 13,  8, 19,  8, 17,  3,  4, 11, 19, 17, 17, 10, 19, 15, 15, 15,  7,  8,  7,  1, 17, 16, 18,  4, 18,  8, 11,  8,  8,  4,  9, 12, 16,  3, 13,  6, 19, 13,  9,  5, 15,  8
    503. , 11, 12, 14,  5,  1, 15, 16,  7,  4, 16,  4,  0,  7,  7,  7,  8, 14,  3, 10,  3,  0,  4,  2, 12, 18,  1, 15,  6,  4,  6,  5, 18,  0,  9, 16,  2, 17, 15,  9,  7, 15,  8, 19,  0, 16, 15, 16,  0,  2, 17, 12,  0,  6,  3, 11,  2, 10, 18, 16, 18, 16, 18,  6,  6
    504. ,  1,  2, 13,  1,  9, 16, 16, 16, 17,  1, 18,  2, 19,  1, 15, 10,  2,  0,  3, 16, 18,  4, 19, 11, 14, 13,  2,  3, 18,  6, 16,  7, 10,  2, 18,  4, 15,  9, 12, 18,  0,  7, 18,  1, 14,  8,  7,  2,  9,  4, 15,  4,  1,  0, 16,  2, 15,  1, 15, 14,  0, 13,  4,  3
    505. , 12, 16, 10, 11,  2,  4, 19, 16,  0,  3,  8,  1,  5, 16,  7, 15, 10,  5,  0, 18,  2, 14, 14,  3, 14, 14,  9, 10,  3,  4,  6, 13,  1,  9, 10, 16,  6, 12,  4,  5, 14,  8,  2, 18, 18,  0,  4, 11,  4,  7,  3, 11, 16,  8,  4,  8, 19,  1,  3,  1, 17, 16, 14,  4
    506. , 19, 14, 12,  9,  6,  0,  9, 17,  5,  0, 15,  9,  1, 10,  4, 10, 15,  2,  4,  4, 14, 12,  6,  8,  0,  4, 18,  5,  5, 10, 10, 13, 13,  8, 13, 18,  0, 12, 13,  0, 12,  1,  5, 14,  7, 11, 11, 12,  1,  3, 10,  4,  4, 10,  6,  0,  7, 12,  4,  2,  9,  1, 18, 16
    507. , 14,  3, 12, 13,  4, 18,  1, 17,  8,  0, 14,  2, 12, 14,  1,  0, 16, 13,  2, 19, 14,  6,  7, 16,  1, 11, 16, 11, 18, 15,  3, 18,  1, 12, 19,  1,  9,  8, 10,  3,  7,  9,  5,  4, 18,  9,  3,  8, 19,  0,  3,  2,  5,  2,  2,  1,  5,  3, 19,  2,  2,  8,  4, 16
    508. , 12, 19, 19,  3,  9,  3,  0,  9, 13, 18, 18, 12,  8,  1,  9, 10,  0,  8, 12, 11,  3,  8, 16, 13,  5,  1, 16, 14,  9,  7, 17,  8,  4, 16, 19,  2, 19, 19, 16, 15,  6,  1,  2,  5,  0,  8, 19, 12,  2,  8,  0,  3,  4, 16, 16,  8,  1, 14, 17,  0, 16, 19, 19,  4
    509. , 12,  1,  8,  2, 17,  9,  9, 18,  6,  7,  9,  5, 13,  1,  5, 14,  7, 19, 18, 14, 14,  0,  1,  5,  9,  0,  2,  5,  8, 14,  7,  6, 19,  0, 15, 16, 10,  7, 18,  9,  9, 11,  5, 18,  7, 13,  6, 12,  3,  4, 18,  7,  5,  1, 12, 11, 19,  3,  6, 18,  2, 15, 19,  7
    510. ,  0,  2, 10,  3,  9, 14,  6,  4,  0,  4,  2, 17,  3,  6, 18, 16,  8,  8,  1, 13, 14,  4, 11,  1,  0, 17, 19,  4, 18, 16, 13, 19, 17, 19, 15, 17,  7,  2,  7,  3,  7,  0,  7,  9, 19,  4,  4,  1,  5,  8, 11,  5, 16,  8, 14,  7, 17,  8, 18, 17,  5,  2, 18,  3
    511. ,  9,  3,  0,  3, 17,  6, 15,  2, 17,  7,  2,  6, 14, 17,  8, 14,  4, 17, 15,  2,  9, 18, 16, 16,  2, 19,  1, 18, 19, 19, 13,  2,  4, 13,  2, 19,  8,  9,  2,  0,  5, 10, 16,  6,  3,  4,  8,  6, 17,  0, 13, 11, 16,  8, 17,  3,  6,  3,  5, 16,  6, 13, 12, 19
    512. , 10, 18,  4, 16,  5,  7,  7,  9,  7, 11, 18,  9, 11,  7,  1,  4, 17,  3,  6,  3, 10,  5, 11, 14,  5,  4, 18, 17,  8, 19, 17,  1,  9,  4, 18,  5,  0, 10,  5,  2,  2, 16,  2, 10,  6,  2, 15, 12, 18, 17, 11, 11,  6, 10, 18,  9, 14, 19, 17, 14,  2,  7, 13, 14
    513. ,  3, 19,  4,  6, 12, 11, 18, 16,  8,  0, 10, 14, 16, 19,  0, 13, 10,  4,  4, 18,  3,  5, 18,  9,  8, 18, 19,  8,  0,  1, 16, 16,  6, 17,  8,  4,  1, 10, 17,  6,  3,  3,  8, 13, 10, 12, 15, 15,  0, 10, 13, 17, 18, 12,  5, 12,  7, 10, 14, 15, 10,  3, 19,  4
    514. ,  8,  9,  1,  9, 16,  1,  6,  1,  5, 14, 17, 14, 10,  8,  0, 13,  3, 11,  6,  6,  4, 10, 15,  7, 14, 16, 19, 19,  1, 19, 10,  6, 16, 16,  5, 13, 17,  8, 18, 19, 16,  4,  7, 17, 13,  5,  3,  0, 16, 17, 12,  0,  5, 13,  2, 18, 15,  1, 19,  7,  4,  3,  7, 10
    515. ,  6, 15, 17, 11,  9, 11, 18, 16,  4,  2, 14,  2, 10, 10, 17,  9, 19, 19,  5, 16,  6, 10,  3, 17,  7, 13, 13, 17, 16, 10, 15, 17,  4,  6, 10, 10, 19,  5, 17,  9, 19, 17, 12, 10,  6, 18,  2, 19,  3,  3, 15, 15, 13, 18,  8, 16, 12, 14,  8, 10, 18,  1,  2, 10
    516. , 16,  0, 14,  0,  5, 17,  4,  4, 14, 15,  7,  4, 13,  0,  1, 19, 15, 17, 18,  7, 13, 13, 18,  8,  6, 19,  2,  1, 16,  6, 17,  6,  6,  8, 17,  5,  9, 11, 15,  7, 12, 12, 10,  3, 11, 11,  3,  4,  3,  5, 12,  0, 10,  8, 17, 14,  2,  4,  2, 13,  7, 16, 15,  9
    517. , 15,  2,  2, 19, 18,  4,  8,  0, 12,  2, 13, 17, 17,  0,  1,  4,  0, 17,  0, 10,  1, 13,  1,  4, 19, 17,  4,  9,  6, 16,  4,  6, 18,  7,  5,  1, 19,  5, 16, 18, 11,  4,  5, 12, 17,  2, 10,  4, 13,  3,  4, 19, 14,  7, 14, 18, 14,  9, 14, 11,  0, 18,  7, 18
    518. ,  5, 18, 14, 11,  7, 15, 12,  5,  2,  2,  1,  4,  9, 10,  3, 15,  7, 10,  9,  2,  9,  8, 12, 16,  0, 19, 13,  4, 17, 16,  6,  8,  7,  6,  1, 18,  8,  1, 17, 15,  0, 15,  9,  3,  1,  4,  7, 19, 13,  4,  2,  9, 11,  1, 13, 13,  4, 19,  5,  8, 16,  0,  5,  9
    519. , 10, 18,  9, 10,  6, 10, 16, 18, 14, 11, 11,  2, 11, 13, 12, 18,  2, 19, 16, 18, 10, 13, 19, 19, 15, 15,  2, 18,  8,  5, 10, 17,  5,  1, 12,  6, 19, 13,  7, 18,  0,  7, 10,  3,  4,  9,  3, 16,  0,  0,  3, 13,  8, 10,  9,  8, 15, 14,  8,  5, 14,  4, 19, 16
    520. ,  4, 15, 14, 13, 12, 13, 15,  6,  1,  9,  6,  6, 19,  5,  1, 12, 14, 15,  2,  4, 16, 18,  1,  2, 16, 17, 19,  5,  4, 13, 10,  5,  1, 18,  6, 16, 17,  2,  3, 15,  0,  6,  5,  5, 15,  9, 17, 10, 12,  2,  8, 18,  1,  3,  5,  4, 18,  7, 16,  4, 14,  4, 10,  3
    521. , 10,  2, 15,  6,  9,  6, 13,  0,  1, 11, 17, 10, 12,  0, 17, 14,  1, 15, 17, 15,  6,  0,  9, 19, 10,  7,  8,  0,  1, 17, 19,  9, 19,  8, 19, 17, 14,  0, 12,  2,  7, 18,  2,  3,  0, 14,  7,  5,  9, 17,  7,  0,  4,  8,  9, 13,  9,  0, 15,  8, 14, 10, 16, 16
    522. , 14, 14,  2,  6, 18, 15,  7,  6,  0, 16,  6, 18,  4,  0, 19,  5,  4, 15, 15,  9, 12, 12,  8, 19,  7,  2,  9, 10, 10,  8,  5, 11,  5,  1, 13,  2,  0,  4, 15,  1,  1, 10,  2, 17, 16, 17,  3, 13,  9, 10,  7, 17, 13,  3, 14, 11,  0,  3, 16, 12,  9, 12,  8,  6
    523. ,  4, 11, 18, 13, 13, 14, 11, 12,  5,  8,  1,  5,  8, 13, 14,  3,  8,  7,  9, 12,  4, 13, 10, 16, 18,  7,  2,  5, 17, 18, 17, 15, 16, 17,  7,  3, 12, 15,  1, 10, 15,  7,  1,  9,  9, 17,  4, 15, 12,  3, 15, 14, 14,  8,  4,  4, 16,  7, 14, 18, 13, 15,  4,  9
    524. , 10, 12, 13,  3,  8, 12,  8, 14, 13, 16, 18,  8, 18, 18, 11,  0, 18,  8,  7, 18,  5,  0,  3, 15,  9,  3,  0,  2,  6, 19,  9,  7, 18, 15, 18, 15,  2,  1, 10, 10,  1, 16, 19, 14,  7, 11, 19,  8, 19,  8, 16,  6, 14, 19, 13,  9, 11, 18,  2, 16,  7,  5,  4,  4
    525. , 13, 16, 16, 12, 12,  5, 12,  4, 10,  9,  9, 10, 10, 13, 10, 13, 14,  7, 15,  0, 14, 12,  7,  6,  9,  7,  5,  2,  3, 16, 19, 12, 11,  0,  0,  0,  7,  1, 15,  1,  3,  9,  1, 18,  3, 12,  7,  3,  3,  6, 11,  0,  8,  0, 19,  2,  0,  6,  8, 15, 12, 11, 19, 15
    526. , 10, 11,  2,  6,  9,  3,  3,  8,  9,  7,  4, 14, 10, 14,  0, 12,  4,  1,  8,  7,  8,  1,  9,  1, 11,  0, 10, 16,  3,  4, 17, 12,  4, 15,  7,  6, 18, 10,  7, 16,  9,  7, 17, 13, 11, 18, 13,  1, 19, 11, 10, 17,  0, 12, 14, 10, 15, 12, 15, 15, 11, 10, 10,  5
    527. ,  4, 10, 15, 18, 10,  6,  1, 19,  6, 17, 15, 18,  0,  1, 17, 14,  5, 17, 19, 18,  2,  5,  5,  2,  5,  7, 16,  2,  8,  7, 12, 10,  5,  9, 10,  5,  2,  2,  1, 19,  1, 17, 16, 19,  2, 18, 15,  1,  4, 15, 14, 18, 15,  3, 16, 16, 17, 18, 18,  7,  3, 11,  7,  4
    528. ,  1,  6,  7,  4,  0,  1, 19,  3,  1,  9, 11, 13,  3, 11,  0,  3, 10, 16,  0,  1, 18, 14,  4,  5, 18,  9,  6, 10, 13, 17, 10,  3, 12,  3,  3,  5,  3, 17,  9, 14, 18, 13, 19, 15,  0, 19, 11, 13,  9, 15,  3, 12, 10,  4, 14,  1,  9, 12,  0, 12, 10,  8, 10, 18
    529. ,  7, 11, 14,  8,  8,  4,  8,  7, 17,  4,  4,  8,  5, 12, 11,  6,  4, 18, 16, 14, 18,  7, 18,  0,  7, 19,  3,  6, 10, 13,  6, 11, 17,  1,  4, 15,  0, 16,  9,  7,  3, 11,  2,  0, 15,  5, 19,  0,  1, 19,  1,  9,  2, 14,  0, 15,  0,  8, 11, 15,  2,  1, 10,  9
    530. ,  3,  5,  4,  9, 12,  7, 13,  5,  8,  3, 11,  7,  0,  8,  1, 19, 17,  4, 15,  8,  0, 11,  9,  8, 13,  4,  4,  2, 12,  5, 18, 11,  2,  4,  9,  9, 14, 17, 17, 11, 12, 18, 18, 19,  5, 14,  5, 18,  4, 15,  8, 10, 12, 14,  6, 18, 10,  5, 17,  6, 18, 19,  0, 16
    531. ,  3,  9,  1, 13, 14,  9,  8,  7,  5,  1, 16,  1,  0, 16, 14,  0, 17, 18, 16,  7,  4, 11,  3, 19,  6,  4,  8, 15, 15,  3,  2,  3, 10,  7,  3, 17,  7,  3,  4, 19,  7, 13, 15, 11, 19,  5,  1,  7, 10, 12, 16,  4, 16, 15, 18, 19, 12, 14, 12,  6,  3,  8, 16, 16
    532. , 11,  5,  0, 14, 17, 16,  6, 11,  6,  5, 18,  1, 19,  8, 19,  6,  6,  8,  0,  2, 11, 12,  8, 12, 12,  1,  6, 12, 15,  0, 19,  4,  4, 19, 16, 10,  5, 13, 15,  8,  3,  1,  1, 13,  0, 18,  7, 14, 10, 15, 18,  9,  3, 12, 13,  1,  7,  0,  3,  6,  5, 10,  9,  5
    533. ,  1,  9,  2, 13,  0, 12,  6, 12, 10,  7,  3,  6,  5,  2, 17, 14, 11, 11,  2,  9,  4,  1, 19,  2,  3,  5, 17, 18,  0, 16,  3,  3, 13, 13,  0, 12,  9,  9, 12, 19,  3, 19,  4,  9,  1,  4, 10, 10,  7,  2, 19,  1,  5,  2, 19, 11, 11, 16, 17,  4, 16,  2,  4, 19
    534. , 15,  1,  8, 16,  5,  2,  9, 15, 14, 19,  9, 14,  5, 17,  0, 13,  4,  8, 17,  4,  7,  3,  0,  8,  4,  8,  0, 17, 10, 17,  3,  5,  3,  4,  0,  2, 17, 10,  3,  8,  6, 11, 15, 15, 19, 15, 12, 15,  2, 19, 13,  1,  9,  1, 13, 13, 19,  6,  0, 10,  9,  9,  9,  1
    535. ,  7, 14,  1, 19, 14, 17, 11, 13,  3,  5, 18, 17, 16,  3,  7,  8,  7,  8, 12, 15,  3, 10,  3,  0, 18, 11, 13, 11, 13, 12, 15, 12,  4,  2,  3,  8,  7,  7, 15, 16, 11, 10, 14,  3,  1,  8, 16, 18, 19, 13,  0, 15,  3, 12,  8,  8, 15,  5, 19, 16, 16,  7, 12,  4
    536. , 19, 17,  1, 14,  8,  8, 11,  8,  5,  5,  2,  4, 19, 13, 11, 12,  3,  4,  0,  4, 11,  4,  2,  3,  7,  5, 11, 11, 17,  0, 15,  0, 19,  9, 13, 18,  0, 17, 14,  6,  0, 17, 18, 12,  9, 10,  4,  9,  1,  1, 15, 17, 18,  7, 15, 11,  8, 13, 18,  0,  3,  8, 18,  7
    537. , 19,  4,  3, 12, 11,  7, 11,  6,  6,  9,  2, 19, 18, 19, 16,  2,  7,  9,  6,  1, 16,  4,  5,  4,  5, 16, 16,  6, 18,  5, 13, 10, 14, 11,  8,  1,  4, 13, 14, 14,  8,  0, 15, 10,  2, 12, 16,  3,  5,  9,  3, 18, 12, 10, 10, 12,  0,  4, 15, 18,  4, 18,  4, 10
    538. , 14,  4,  6, 11,  9,  5,  7, 15,  0, 17,  1, 12,  4,  1,  5,  7,  9, 12,  3,  0,  8, 10,  2, 16,  1,  8,  8, 10, 12, 13, 18,  8,  7,  1, 10,  3,  8,  3,  8, 19,  0, 12,  3,  4, 14, 14, 15, 12,  2,  1, 12,  7, 10,  0,  9,  6,  9,  2,  8,  4, 18, 17,  8, 10
    539. ,  0, 10, 17, 18,  0, 17, 16, 11,  7, 14,  4,  4,  2,  6, 10, 13,  6, 16, 11, 16,  4,  6,  2, 16, 12, 14, 17, 18,  5,  2,  8, 17, 14, 13,  9,  5,  9, 14,  4, 13, 19, 14, 16, 14,  0,  6, 18, 13, 19, 13,  8, 15, 10,  9, 16,  4,  9, 12,  5,  0,  8, 19, 17,  1
    540. , 13,  0,  9, 11,  4,  6,  2, 18,  4, 10, 17, 15, 17,  0,  4, 16, 15,  3,  2,  2,  8,  0,  1,  8, 11,  7,  3,  9, 12, 18, 16, 14, 18, 13,  8,  4, 13, 11,  4,  9,  2, 10, 16,  1, 15, 18, 19,  1, 11, 13,  8, 11, 12,  6,  4, 11, 19, 17, 18,  8,  1,  0,  0, 18
    541. ,  4, 17, 15,  9,  2,  3, 19, 11,  8, 15,  4, 15,  8,  7,  2,  0,  6, 13, 10, 15, 19,  7,  5,  1, 19, 17,  6, 14,  7, 15, 12,  2, 14,  4, 15, 18,  9,  0, 16, 11,  0, 15, 17,  9,  0, 10, 12,  7, 11, 19, 15,  8,  0,  9,  9, 19, 12,  9,  8,  8, 18, 11,  7, 19
    542. , 19,  7, 12,  5,  5, 14, 14,  7, 17, 19, 16, 16,  1, 12,  1,  0, 10,  6, 18,  1,  1, 12,  3, 14,  3,  8,  3, 19, 10,  1, 14,  4,  9, 19, 14,  7,  0,  3,  7, 18,  6, 12, 18, 12,  8,  5, 14,  0, 16,  6,  5, 13,  4,  2, 12, 17,  9,  8,  7,  9,  6, 11,  7,  6
    543. , 19,  1,  4,  4,  0,  1, 16,  6, 19,  5, 14,  3,  3,  8, 10,  5, 19, 19, 16, 15,  3,  4, 19, 17,  6, 18,  5, 17, 14, 19,  8,  2, 14,  5,  8, 16, 15, 16, 14,  2,  8, 15, 18,  0, 11, 17, 12,  3, 17,  0, 19, 18, 15,  8,  5, 18,  8,  7,  7, 19,  2,  8, 16,  2
    544. ,  7, 18,  4,  9,  5, 16, 19,  3, 10,  9,  5, 11, 17,  6, 14, 16, 19, 13, 18, 14,  1,  2,  2,  0, 18, 17, 16, 14, 15,  7, 10, 13, 11,  8,  5,  4,  8, 12, 18, 16, 15, 15,  7, 12, 15,  6,  6,  6,  4, 10, 16,  0, 18,  4,  0,  8,  8,  9, 19,  9,  7,  5, 15, 10
    545. ,  6, 16, 13,  3,  1,  4,  5, 17,  6, 15,  6,  9, 10, 15, 15,  2, 19,  9, 16,  0, 17,  9,  2, 16,  2,  5,  6,  2, 10,  4, 18,  7,  0, 16, 14, 14, 14,  9, 13,  7, 12, 11,  3, 10,  2, 18,  3,  5, 16,  9, 16,  3,  4, 18,  8,  1, 18,  6,  2,  7, 15, 11, 11,  2
    546. ,  3,  5,  8,  4,  0, 13,  7, 16,  2,  0, 19,  3,  7,  7,  3,  8,  7,  5, 18, 13, 16,  1,  8, 19, 15,  2, 13,  7,  3,  3,  1, 16, 18,  0,  4,  4, 10, 12, 15,  5, 11, 10, 11,  8,  1, 19,  8, 10,  2,  9,  7,  8, 18, 17, 19,  0, 11, 11,  8,  5, 11,  2, 17,  1
    547. ,  0,  3, 14,  4,  7, 10,  6, 17,  7, 10, 15, 10,  0,  9,  7,  9,  6, 15,  6,  4, 14, 18,  4, 19, 19, 18, 12, 13, 19,  7,  2, 15,  7,  5, 19, 10, 16,  8,  4,  4, 19, 10,  7, 10, 10,  0, 16,  9,  4,  9, 12, 18,  4,  8, 17,  0,  7,  7, 16, 15, 11, 17, 18, 18
    548. , 12, 14,  1,  1,  8, 10,  8,  4, 12, 19, 18,  1,  2, 19, 12, 11,  5,  8,  6,  3,  4, 16, 16,  4,  7,  3, 19, 14,  4, 10, 10,  9, 18,  9, 16,  3, 16,  6,  9,  4, 15,  5,  4, 18,  9, 16, 16,  5, 19,  1,  4,  7, 10, 10,  1, 18, 19,  6,  2, 10,  2,  1, 18, 14
    549.     };
    550.     */
    551.     /**
    552.     float A[441] = {
    553.                     13, 20, 16,  7, 29,  0, 23,  4,  6,  2,  1, 15, 21,  3, 27,  1, 20, 23,  6, 21,  4
    554.                     , 20,  1,  8, 16, 24, 28, 18,  8,  4, 26, 13, 19,  6, 12, 24, 27,  6, 22,  6, 17,  9
    555.                     , 16,  8, 21, 11, 20,  9, 11, 10, 13,  7, 22, 10,  7, 19,  8, 28, 26, 20, 22, 19, 16
    556.                     ,  7, 16, 11, 25, 16, 13, 26,  2, 13, 18,  2, 18,  9, 14,  1, 27, 19,  6,  0, 18, 26
    557.                     , 29, 24, 20, 16, 23,  6, 26,  1, 10, 27,  6, 16, 14, 20, 23, 29, 24,  0, 17,  9,  1
    558.                     ,  0, 28,  9, 13,  6, 18,  4, 27, 18, 28, 20, 13,  2, 20, 10,  4,  6,  6,  5, 11, 25
    559.                     , 23, 18, 11, 26, 26,  4, 16,  8,  7, 28, 26, 10, 19, 24, 15,  2, 23, 23, 24, 22, 29
    560.                     ,  4,  8, 10,  2,  1, 27,  8, 28, 20, 17, 18,  1, 26, 25, 25, 16, 22, 15, 15,  1, 24
    561.                     ,  6,  4, 13, 13, 10, 18,  7, 20, 15, 27, 23,  2, 23, 15, 14, 26, 20, 14,  9, 14, 15
    562.                     ,  2, 26,  7, 18, 27, 28, 28, 17, 27,  7,  5, 11, 28, 18, 25, 23,  5,  6, 12, 12, 12
    563.                     ,  1, 13, 22,  2,  6, 20, 26, 18, 23,  5,  7, 22, 13, 14, 11,  0, 17, 23,  3,  8, 29
    564.                     , 15, 19, 10, 18, 16, 13, 10,  1,  2, 11, 22, 25, 13, 25, 12, 18, 14,  7, 23,  2, 16
    565.                     , 21,  6,  7,  9, 14,  2, 19, 26, 23, 28, 13, 13,  3, 10, 25,  6, 26, 28, 14,  9,  5
    566.                     ,  3, 12, 19, 14, 20, 20, 24, 25, 15, 18, 14, 25, 10,  0,  4, 29, 23,  0, 24,  6, 23
    567.                     , 27, 24,  8,  1, 23, 10, 15, 25, 14, 25, 11, 12, 25,  4, 10,  9,  9,  7,  5, 22,  2
    568.                     ,  1, 27, 28, 27, 29,  4,  2, 16, 26, 23,  0, 18,  6, 29,  9, 14, 22, 19,  8, 16, 29
    569.                     , 20,  6, 26, 19, 24,  6, 23, 22, 20,  5, 17, 14, 26, 23,  9, 22, 28, 20, 25, 29, 13
    570.                     , 23, 22, 20,  6,  0,  6, 23, 15, 14,  6, 23,  7, 28,  0,  7, 19, 20,  5, 15, 13, 13
    571.                     ,  6,  6, 22,  0, 17,  5, 24, 15,  9, 12,  3, 23, 14, 24,  5,  8, 25, 15, 17, 27, 18
    572.                     , 21, 17, 19, 18,  9, 11, 22,  1, 14, 12,  8,  2,  9,  6, 22, 16, 29, 13, 27,  3, 23
    573.                     ,  4,  9, 16, 26,  1, 25, 29, 24, 15, 12, 29, 16,  5, 23,  2, 29, 13, 13, 18, 23, 22
    574.     };
    575.     */
    576.     printf("\nA=\n");
    577.     printFloatMatrixRowMajor(A, m, n, n);
    578.  
    579.     int jt;
    580.     float eps;
    581.     int result = 0;
    582.     float v_A[49];
    583.  
    584.     jt = 89;
    585.     eps = 0.000001;
    586.  
    587.     //result = jcbi(A, n, v_A, eps, jt);
    588.     result = jcbi_roll(A, n, v_A, eps, jt);
    589.     
    590.     printf("\nSigma_A=\n");
    591.     printFloatMatrixRowMajor(A, n, n, n);
    592.     printf("\nV_A=\n");
    593.     printFloatMatrixRowMajor(v_A, n, n, n);
    594.  
    595. }
    596.  
    597.  
    598. int main() {
    599.     svd_sym_Test();
    600.     //svd_Little_mGEn_Test();// A'*A => V S
    602.     return 0;
    603. }

  • 相关阅读:
    EMAS Serverless 到底有多便利？
    自动驾驶学习笔记（八）——路线规划
    图像滤波算法
    Vuex环境搭建
    全球第4大操作系统(鸿蒙)的软件后缀.hap
    大话机器学习准确率(Accuracy)、精确率(Pecision)、召回率(Recall)以及TP、FP、TN、FN
    H3C 6520X 配置IRF+BFD
    C专家编程 第7章 对内存的思考 7.4 cache存储器
    【深度学习笔记】计算机视觉——图像增广
    Linux基础操作
  • 原文地址：https://blog.csdn.net/eloudy/article/details/128116759