xunhuan Jacobi method SVD

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

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

void printFloatMatrixRowMajor(float* A, int M, int N, int lda);
void printFloatMatrixColMajor(float* A, int M, int N, int lda);



int jcbi_roll(float* a, int n, float* v, float eps, int jt) {
        int i, j, p, q, u, w, t, s, loops;
        float fm, cn, sn, omega, x, y, d;
        loops = 1;
        int contiCount = 0;

        for (i = 0; i <= n - 1; i++)//Let V = I;
        {
            v[i * n + i] = 1.0;
            for (j = 0; j <= n - 1; j++)
                if (i != j) v[i * n + j] = 0.0;
        }

        while (1 == 1)
        {
            fm = 0.0;
            for (p = 1; p <= n - 1; p++)
            {
                for (q = 0; q <= p - 1; q++)
                {
                    if (loops > jt) { printf("\n\ncontiCount = %d\n", contiCount);  return(-1); }// 超过迭代次数,依然没有收敛,退出,返回-1

                    loops = loops + 1;//loopCount++ ;
                    u = p * n + q; //A(p, q)
                    w = p * n + p; //A(p, p)
                    t = q * n + p; //A(q, p)
                    s = q * n + q; //A(q, q)
                    if (fabs(a[u]) < eps) { contiCount++; continue; }

                    //计算sin(2θ)
                    x = -a[u];                          // m = -A(p,q)
                    y = (a[s] - a[w]) / 2.0;            // n = 1/2(A(q,q)-A(pp))
                    omega = x / sqrt(x * x + y * y);    // m/(sqrt(m^2 + n^2))             
                    if (y < 0.0) omega = -omega;        //sgn(n)*              sin(2θ)=omega
                    //计算sin(θ)
                    sn = 1.0 + sqrt(1.0 - omega * omega);
                    sn = omega / sqrt(2.0 * sn);        // sin(θ)
                    //计算cos(θ)
                    cn = sqrt(1.0 - sn * sn);           // cos(θ)
                    fm = a[w];                          // A(p,p)
                    //计算 四个元素
                    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θ)
                    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θ)
                    a[u] = 0.0;                                          // A(p,q)' = 0;
                    a[t] = 0.0;                                          // A(q,p)' = 0;
                    //计算第p行 和第q行
                    for (j = 0; j <= n - 1; j++)            // updat row p  and  row q;
                        if ((j != p) && (j != q))
                        {
                            u = p * n + j;               // 遍历第 p 行的每一个元素
                            w = q * n + j;               // 遍历第 q 行的每一个元素
                            fm = a[u];                    // 
                            a[u] = fm * cn + a[w] * sn;    // 自己的上次值 累加上 伴随行的元素上次值
                            a[w] = -fm * sn + a[w] * cn;    // 自己的原值 累加上 伴随行的元素值
                        }
                    //计算第p列和第q列
                    for (i = 0; i <= n - 1; i++)            // update column p  and  column q
                        if ((i != p) && (i != q))
                        {
                            u = i * n + p;               // A(i, p) 遍历 p-th column 
                            w = i * n + q;               // A(i, q) 遍历 q-th column
                            fm = a[u];                    //
                            a[u] = fm * cn + a[w] * sn;     // 自己的上次值 累加上 伴随行的元素上次值
                            a[w] = -fm * sn + a[w] * cn;    // 自己的上次值 累加上 伴随行的元素上次值
                        }
                    for (i = 0; i <= n - 1; i++)            // update EigenVectors p-th  and q-th
                    {
                        u = i * n + p;                      // V(i, p) 遍历 p-th column
                        w = i * n + q;                      // V(i, q) 遍历 q-th column
                        fm = v[u];                          //
                        v[u] = fm * cn + v[w] * sn;         // 
                        v[w] = -fm * sn + v[w] * cn;        // 
                    } 
                }///LLLL
            }
        }

        printf("\n\ncontiCount=%d\n", contiCount);

        return(1);
    }


int jcbi(float* a, int n, float* v, float eps, int jt){
    int i, j, p, q, u, w, t, s, loops;
    float fm, cn, sn, omega, x, y, d;
    loops = 1;

    for (i = 0; i <= n - 1; i++)//Let V = I;
    {
        v[i * n + i] = 1.0;
        for (j = 0; j <= n - 1; j++)
            if (i != j) v[i * n + j] = 0.0;
    }

    while (1 == 1)
    {
        fm = 0.0;
        for (i = 1; i <= n - 1; i++) {//only check the lower triangle matrix of D_i, without the diag elements
            for (j = 0; j <= i - 1; j++)// q 恒小于 p;   q 不可能等于 p
            {
                d = fabs(a[i * n + j]);
                if ((i != j) && (d > fm))
                {
                    fm = d;//当前最大值
                    p  = i;//行号
                    q  = j;//列号
                }
            }
        }

        if (fm < eps)  return(1);//the max element is less than eps, tolerance
        if ( loops > jt)  return(-1);// 超过迭代次数,依然没有收敛,退出,返回-1

        loops = loops + 1;//loopCount++ ;
        u = p * n + q; //A(p, q)
        w = p * n + p; //A(p, p)
        t = q * n + p; //A(q, p)
        s = q * n + q; //A(q, q)

//计算sin(2θ)
        x = -a[u];                          // m = -A(p,q)
        y = (a[s] - a[w]) / 2.0;            // n = 1/2(A(q,q)-A(pp))
        omega = x / sqrt(x * x + y * y);    // m/(sqrt(m^2 + n^2))             
        if (y < 0.0) omega = -omega;        //sgn(n)*              sin(2θ)=omega
//计算sin(θ)
        sn = 1.0 + sqrt(1.0 - omega * omega);
        sn = omega / sqrt(2.0 * sn);        // sin(θ)
//计算cos(θ)
        cn = sqrt(1.0 - sn * sn);           // cos(θ)
        fm = a[w];                          // A(p,p)
//计算 四个元素
        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θ)
        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θ)
        a[u] = 0.0;                                          // A(p,q)' = 0;
        a[t] = 0.0;                                          // A(q,p)' = 0;
//计算第p行 和第q行
        for (j = 0; j <= n - 1; j++)            // updat row p  and  row q;
            if ((j != p) && (j != q))
            {
                u    = p * n + j;               // 遍历第 p 行的每一个元素
                w    = q * n + j;               // 遍历第 q 行的每一个元素
                fm   = a[u];                    // 
                a[u] =  fm * cn + a[w] * sn;    // 自己的上次值 累加上 伴随行的元素上次值
                a[w] = -fm * sn + a[w] * cn;    // 自己的原值 累加上 伴随行的元素值
            }
//计算第p列和第q列
        for (i = 0; i <= n - 1; i++)            // update column p  and  column q
            if ((i != p) && (i != q))
            {
                u    = i * n + p;               // A(i, p) 遍历 p-th column 
                w    = i * n + q;               // A(i, q) 遍历 q-th column
                fm   = a[u];                    //
                a[u] = fm * cn + a[w] * sn;     // 自己的上次值 累加上 伴随行的元素上次值
                a[w] = -fm * sn + a[w] * cn;    // 自己的上次值 累加上 伴随行的元素上次值
            }
        for (i = 0; i <= n - 1; i++)            // update EigenVectors p-th  and q-th
        {
            u = i * n + p;                      // V(i, p) 遍历 p-th column
            w = i * n + q;                      // V(i, q) 遍历 q-th column
            fm = v[u];                          //
            v[u] = fm * cn + v[w] * sn;         // 
            v[w] = -fm * sn + v[w] * cn;        // 
        }
    }

    return(1);
}

void transpositionColMajor(float* A, int M, int N, float* At, int lda, int ldat) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < M; j++) {
            At[i + j * ldat] = A[j + i * lda];
        }
    }
}
void transpositionRowMajor(float* A, int M, int N, float* At, int lda, int ldat) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < M; j++) {
            At[i * ldat + j] = A[j * lda + i];
        }
    }
}
void SGemmColMajor(float* A, float* B, int M, int N, int K, float* C, int lda, int ldb, int ldc) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            C[i + j * ldc] = 0.0f;
            for (int k = 0; k < K; k++) {
                C[i + j * ldc] += A[i + k * lda] * B[k + j * ldb];
            }
        }
    }
}
void SGemmRowMajor(float* A, float* B, int M, int N, int K, float* C, int lda, int ldb, int ldc) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            C[i * ldc + j] = 0.0f;
            for (int k = 0; k < K; k++) {
                C[i * ldc + j] += A[i * lda + k] * B[k * ldb + j];
            }
        }
    }
}
void printFloatMatrixColMajor(float* A, int M, int N, int lda) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            printf(" %9.7f", A[i + j * lda]);
        }
        printf("\n\n");
    }
}
void printFloatMatrixRowMajor(float* A, int M, int N, int lda) {
    for (int i = 0; i < M; i++) {
        for (int j = 0; j < N; j++) {
            printf(" %6.4f", A[i * lda + j]);
        }
        printf("\n\n");
    }
}
void gemmTest() {
    int m, n, k;
    m = 4;
    n = 2;
    k = 3;
    float A[12] = { 5,3,6,0,2,8,0,1,4,2,4,6 };
    float B[6] = { 2,1,3,4,3,2 };
    float C[8] = { 0.0 };
    printf("\nA=\n");
    printFloatMatrixColMajor(A, m, k, m);
    printf("\nB=\n");
    printFloatMatrixColMajor(B, k, n, k);
    SGemmColMajor(A, B, m, n, k, C, m, k, m);
    printf("\nC=\n");
    printFloatMatrixColMajor(C, m, n, m);
}
void transpositionTest() {
    int m, n, k;
    m = 4;
    n = 3;

    float  A[12] = { 5,3,6,0,2,8,0,1,4,2,4,6 };
    float At[12];
    //void transpositionColMajor(float* A, int M, int N, float* At, int lda, int ldat)
    transpositionColMajor(A, m, n, At, m, n);

    printf("\nA=\n");
    printFloatMatrixColMajor(A, m, n, m);
    printf("\nAt=\n");
    printFloatMatrixColMajor(At, n, m, n);
}
void svd_sym_Little_Test() {
    int n = 6;

    //float  A[36] = { 2, -1, 0, -1, 2, -1, 0, -1, 2 };
    float A[36] = {
   1.9757,   2.4581,   1.4726, 1.3031, 1.5231, 0.8238,
   2.4581,   3.5969,   1.9414, 2.3095, 2.3901, 1.4690,
   1.4726,   1.9414,   1.5127, 1.2891, 0.9270, 0.5674,
   1.3031,   2.3095,   1.2891, 1.8225, 1.4753, 1.0712,
   1.5231,   2.3901,   0.9270, 1.4753, 2.0327, 0.9849,
   0.8238,   1.4690,   0.5674, 1.0712, 0.9849, 0.9280
    };


    printf("\nA=\n");
    printFloatMatrixRowMajor(A, n, n, n);

    int   jt;
    float eps;
    int result = 0;
    float eigenVs[36];

    jt = 100;
    eps = 0.000001;
    result = jcbi(A, n, eigenVs, eps, jt);

    printf("\neigenValus=\n");
    printFloatMatrixRowMajor(A, n, n, n);


    printf("\neigenVs=\n");
    printFloatMatrixRowMajor(eigenVs, n, n, n);
}
void svd_Little_mGEn_Test() {
    int m, n, k;
    m = 4;
    n = 3;

    //float  A[12] = { 5,3,6,0,2,8,0,1,4,2,4,6 };
    float  A[12] = { 5,2,4,3,8,2,6,0,4,0,1,6 };

    float At[12];
    //void transpositionRowMajor(float* A, int M, int N, float* At, int lda, int ldat)
    transpositionRowMajor(A, m, n, At, n, m);

    printf("\nA=\n");
    printFloatMatrixRowMajor(A, m, n, n);
    printf("\nAt=\n");
    printFloatMatrixRowMajor(At, n, m, m);

    float A_x_At[16];
    float At_x_A[9];

    //void SGemmColMajor(float* A, float* B, int M, int N, int K, float* C, int lda, int ldb, int ldc)
    SGemmRowMajor(At, A, n, n, m, At_x_A, m, n, n);
    printf("\nAt_x_A=\n");
    printFloatMatrixRowMajor(At_x_A, n, n, n);

    int jt;
    float eps;
    int result = 0;
    float v_At_x_A[9];

    jt = 100;
    eps = 0.000001;

    result = jcbi(At_x_A, n, v_At_x_A, eps, jt);
    printf("\nSigma in At_x_A=\n");
    printFloatMatrixRowMajor(At_x_A, n, n, n);
    printf("\nV of At_x_A=\n");
    printFloatMatrixRowMajor(v_At_x_A, n, n, n);

}
void svd_Little_mLTn_Test() {
    int m, n, k;
    m = 4;
    n = 3;

    //float  A[12] = { 5,3,6,0,2,8,0,1,4,2,4,6 };
    float  A[12] = { 5,2,4,3,8,2,6,0,4,0,1,6 };

    float At[12];
    //void transpositionColMajor(float* A, int M, int N, float* At, int lda, int ldat)
    transpositionColMajor(A, m, n, At, m, n);

    printf("\nA=\n");
    printFloatMatrixColMajor(A, m, n, m);
    printf("\nAt=\n");
    printFloatMatrixColMajor(At, n, m, n);

    float A_x_At[16];
    float At_x_A[9];

    //void SGemmColMajor(float* A, float* B, int M, int N, int K, float* C, int lda, int ldb, int ldc)
    /****
    SGemmColMajor(A, At, m, m, n, A_x_At, m, n, m);
    printf("\nA_x_At=\n");
    printFloatMatrixColMajor(A_x_At, m, m, m);
    *****/
    SGemmColMajor(At, A, n, n, m, At_x_A, n, m, n);
    printf("\nAt_x_A=\n");
    printFloatMatrixColMajor(At_x_A, n, n, n);

    int jt;
    float eps;
    int result = 0;
    float v_A_x_At[16];
    float v_At_x_A[9];

    jt = 100;
    eps = 0.000001;

    /****
    result = jcbi(A_x_At, m, v_A_x_At, eps, jt);

    printf("\nA_x_At=\n");
    printFloatMatrixColMajor(A_x_At, m, m, m);
    printf("\nV of A_x_At=\n");
    printFloatMatrixColMajor(v_A_x_At, m, m, m);
    *****/

    result = jcbi(At_x_A, n, v_At_x_A, eps, jt);
    printf("\nAt_x_A=\n");
    printFloatMatrixColMajor(At_x_A, n, n, n);
    printf("\nV of At_x_A=\n");
    printFloatMatrixColMajor(v_At_x_A, n, n, n);

}
void printFloatVector(float* X, int n) {
    for (int idx = 0; idx < n; idx++) {
        printf(" %9.7e,", X[idx]);
    }
    printf("\n");
}
void svd_Big_Test() {
    int m, n, k;
    int batchSize = 1;
    m = 300;
    n = 181;

    float* A;
    float* At;
    A = (float*)malloc(sizeof(float) * m * n * batchSize);
    At = (float*)malloc(sizeof(float) * n * m * batchSize);
    const char* fileName = "test_svd.bin";
    //void transpositionColMajor(float* A, int M, int N, float* At, int lda, int ldat)
    //transpositionColMajor(A, m, n, At, m, n);
    FILE* fp = fopen(fileName, "rb");

    if (fp == NULL) {
        printf("FILE OPEN FAIL!\n");
        return;
    }

    fread(A, sizeof(float), m * n * batchSize, fp);
    //fread(A, sizeof(float), 7, fp);

    transpositionColMajor(A, m, n, At, m, n);

    printf("A=\n\n");
    printFloatVector(A, 7);


    printf("\nA=\n");
    //printFloatMatrixColMajor(A, m, n, m);
    printf("\nAt=\n");
    //printFloatMatrixColMajor(At, n, m, n);

    float* A_x_At;
    float* At_x_A;

    A_x_At = (float*)malloc(m * m * sizeof(float) * batchSize);
    At_x_A = (float*)malloc(n * n * sizeof(float) * batchSize);

    //void SGemmColMajor(float* A, float* B, int M, int N, int K, float* C, int lda, int ldb, int ldc)
    SGemmColMajor(A, At, m, m, n, A_x_At, m, n, m);
    printf("\nA_x_At=\n");
    printFloatMatrixColMajor(A_x_At, m, m, m);

    SGemmColMajor(At, A, n, n, m, At_x_A, n, m, n);
    printf("\nAt_x_A=\n");
    printFloatMatrixColMajor(At_x_A, n, n, n);

    int jt;
    float eps;
    int result = 0;
    float* v_A_x_At;
    float* v_At_x_A;
    v_A_x_At = (float*)malloc(m * m * sizeof(float) * batchSize);
    v_At_x_A = (float*)malloc(n * n * sizeof(float) * batchSize);

    jt = 100;
    eps = 0.000001;
    result = jcbi(A_x_At, m, v_A_x_At, eps, jt);

    printf("\nA_x_At=\n");
    printFloatMatrixColMajor(A_x_At, m, m, m);
    printf("\nU of A_x_At=\n");
    printFloatMatrixColMajor(v_A_x_At, m, m, m);


    result = jcbi(At_x_A, n, v_At_x_A, eps, jt);
    printf("\nAt_x_A=\n");
    printFloatMatrixColMajor(At_x_A, n, n, n);
    printf("\nV of At_x_A=\n");
    printFloatMatrixColMajor(v_At_x_A, n, n, n);
    /****/
}
void svd_sym_Test() {
    int m, n;
    m = 7;
    n = 7;
    float A[49] = {
                                  6, 6, 0, 9, 7, 2, 6,
                                  6, 5, 9, 4, 2, 5, 5,
                                  0, 9, 8, 3, 5, 0, 6,
                                  9, 4, 3, 8, 9, 9, 4,
                                  7, 2, 5, 9, 0, 6, 0,
                                  2, 5, 0, 9, 6, 4, 4,
                                  6, 5, 6, 4, 0, 4, 9
    };


    /**
    float A[4096] = {
   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
, 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
,  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
,  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
,  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
,  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
, 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
, 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
, 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
,  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
,  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
,  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
, 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
,  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
,  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
,  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
, 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
, 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
, 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
,  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
, 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
, 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
, 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
, 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
, 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
,  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
,  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
, 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
,  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
,  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
,  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
, 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
, 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
,  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
, 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
,  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
, 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
, 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
,  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
, 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
, 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
, 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
,  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
,  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
,  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
,  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
,  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
, 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
,  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
, 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
,  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
, 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
, 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
, 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
,  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
, 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
,  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
, 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
, 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
,  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
,  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
,  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
,  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
, 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
    };
    */
    /**
    float A[441] = {
                    13, 20, 16,  7, 29,  0, 23,  4,  6,  2,  1, 15, 21,  3, 27,  1, 20, 23,  6, 21,  4
                    , 20,  1,  8, 16, 24, 28, 18,  8,  4, 26, 13, 19,  6, 12, 24, 27,  6, 22,  6, 17,  9
                    , 16,  8, 21, 11, 20,  9, 11, 10, 13,  7, 22, 10,  7, 19,  8, 28, 26, 20, 22, 19, 16
                    ,  7, 16, 11, 25, 16, 13, 26,  2, 13, 18,  2, 18,  9, 14,  1, 27, 19,  6,  0, 18, 26
                    , 29, 24, 20, 16, 23,  6, 26,  1, 10, 27,  6, 16, 14, 20, 23, 29, 24,  0, 17,  9,  1
                    ,  0, 28,  9, 13,  6, 18,  4, 27, 18, 28, 20, 13,  2, 20, 10,  4,  6,  6,  5, 11, 25
                    , 23, 18, 11, 26, 26,  4, 16,  8,  7, 28, 26, 10, 19, 24, 15,  2, 23, 23, 24, 22, 29
                    ,  4,  8, 10,  2,  1, 27,  8, 28, 20, 17, 18,  1, 26, 25, 25, 16, 22, 15, 15,  1, 24
                    ,  6,  4, 13, 13, 10, 18,  7, 20, 15, 27, 23,  2, 23, 15, 14, 26, 20, 14,  9, 14, 15
                    ,  2, 26,  7, 18, 27, 28, 28, 17, 27,  7,  5, 11, 28, 18, 25, 23,  5,  6, 12, 12, 12
                    ,  1, 13, 22,  2,  6, 20, 26, 18, 23,  5,  7, 22, 13, 14, 11,  0, 17, 23,  3,  8, 29
                    , 15, 19, 10, 18, 16, 13, 10,  1,  2, 11, 22, 25, 13, 25, 12, 18, 14,  7, 23,  2, 16
                    , 21,  6,  7,  9, 14,  2, 19, 26, 23, 28, 13, 13,  3, 10, 25,  6, 26, 28, 14,  9,  5
                    ,  3, 12, 19, 14, 20, 20, 24, 25, 15, 18, 14, 25, 10,  0,  4, 29, 23,  0, 24,  6, 23
                    , 27, 24,  8,  1, 23, 10, 15, 25, 14, 25, 11, 12, 25,  4, 10,  9,  9,  7,  5, 22,  2
                    ,  1, 27, 28, 27, 29,  4,  2, 16, 26, 23,  0, 18,  6, 29,  9, 14, 22, 19,  8, 16, 29
                    , 20,  6, 26, 19, 24,  6, 23, 22, 20,  5, 17, 14, 26, 23,  9, 22, 28, 20, 25, 29, 13
                    , 23, 22, 20,  6,  0,  6, 23, 15, 14,  6, 23,  7, 28,  0,  7, 19, 20,  5, 15, 13, 13
                    ,  6,  6, 22,  0, 17,  5, 24, 15,  9, 12,  3, 23, 14, 24,  5,  8, 25, 15, 17, 27, 18
                    , 21, 17, 19, 18,  9, 11, 22,  1, 14, 12,  8,  2,  9,  6, 22, 16, 29, 13, 27,  3, 23
                    ,  4,  9, 16, 26,  1, 25, 29, 24, 15, 12, 29, 16,  5, 23,  2, 29, 13, 13, 18, 23, 22
    };
    */
    printf("\nA=\n");
    printFloatMatrixRowMajor(A, m, n, n);

    int jt;
    float eps;
    int result = 0;
    float v_A[49];

    jt = 89;
    eps = 0.000001;

    //result = jcbi(A, n, v_A, eps, jt);
    result = jcbi_roll(A, n, v_A, eps, jt);
    
    printf("\nSigma_A=\n");
    printFloatMatrixRowMajor(A, n, n, n);
    printf("\nV_A=\n");
    printFloatMatrixRowMajor(v_A, n, n, n);

}


int main() {
    svd_sym_Test();
    //svd_Little_mGEn_Test();// A'*A => V S

    return 0;
}

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值