マトリックス求逆アルゴリズム-全選択メタガウス-ヨルダン法

マトリックス求逆アルゴリズム-全選択メタガウス-ヨルダン法
Tags:逆行列
全選主元ガウス-ヨルダン法の逆を求める手順は以下の通りである.
1.kについて0からn-1まで以下のステップを行う.

は、k行目、k列目から始まる右下のサブアレイから絶対値が最も大きい要素を選択し、行交換と列交換によってプライマリ要素の位置に交換される行番号と列番号を記憶します.このステップを全選択マスターと呼びます.

m(k, k) = 1/m(k, k)

m(k, j) = m(k, j) * m(k, k),j = 0, 1, ..., n-1;j != k

m(i, j) = m(i, j) - m(i, k) * m(k, j),i, j = 0, 1, ..., n-1;i, j != k

m(i, k) = -m(i, k) * m(k, k),i = 0, 1, ..., n-1;i != k

2.全選択マスタープロセスで記録された行、列によって交換された情報に基づいて回復する.回復の原則は以下の通りである.
すべてのマスターを選択する過程で、先に交換した行(列)の後に回復します.元の行(列)交換は列(行)交換で復元される.
サンプルコード
1.マトリックス構造定義

#ifdef DEBUG
#define LM_ASSERT(x)    do{ if(!(x)) while(1); }while(0)
#else
#define LM_ASSERT(x)
#endif

typedef float real_t;

typedef struct matrix
{
    unsigned int row;
    unsigned int col;
    real_t *m;
}matrix_t;

#define MATRIX(M, r, c) (*((M)->m + (r)*(M)->col + (c)))

2関数実装

#define FABS(x)     fabs(x)

/**
 * @brief    
 */
void matrix_swap_row(matrix_t *m, unsigned int i, unsigned int j)
{
    unsigned int k;
    real_t tmp;

    LM_ASSERT(i < m->row);
    LM_ASSERT(j < m->row);

    for(k=0; kcol; k++)
    {
        tmp = MATRIX(m, i, k);
        MATRIX(m, i, k) = MATRIX(m, j, k);
        MATRIX(m, j, k) = tmp;
    }
}

/**
 * @brief    
 */
void matrix_swap_col(matrix_t *m, unsigned int i, unsigned int j)
{
    unsigned int k;
    real_t tmp;

    LM_ASSERT(i < m->col);
    LM_ASSERT(j < m->col);

    for(k=0; krow; k++)
    {
        tmp = MATRIX(m, k, i);
        MATRIX(m, k, i) = MATRIX(m, k, j);
        MATRIX(m, k, j) = tmp;
    }
}

/**
 * @brief     
 */
void matrix_copy(matrix_t *to, matrix_t *from)
{
    unsigned int i, j;

    matrix_reshape(to, from->row, from->col);

    for(i=0; irow; i++)
        for(j=0; jcol; j++)
            MATRIX(to, i, j) = MATRIX(from, i, j);
}

/**
 * @brief       
 */
int matrix_reshape(matrix_t *m, unsigned int row, unsigned int col)
{
    LM_ASSERT(m != NULL);
    LM_ASSERT(row != 0);
    LM_ASSERT(col != 0);

    if (row * col == m->row * m->col)
    {
        m->row = row;
        m->col = col;
    }
    else
    {
        if (m->m != NULL)
            free(m->m);

        m->m = malloc(row * col * sizeof(real_t));
        if (m->m != NULL)
        {
            m->row = row;
            m->col = col;
        }
        else
        {
            m->row = m->col = 0;
            return -1;
        }
    }

    return 0;
}

/**
 * @brief     
 */
int matrix_inv(matrix_t *inv, matrix_t *a)
{
    int i, j, k;
    int ret = 0;

    //!      
    LM_ASSERT(a->row == a->col);

    unsigned int is[a->row];
    unsigned int js[a->col];

    real_t max;

    matrix_reshape(inv, a->row, a->col);

    matrix_copy(inv, a);

    for(k=0; krow; k++)
    {
        //step 1,     
        max = 0;
        is[k] = k;
        js[k] = k;

        for(i=k; irow; i++)
        {
            for(j=k; jcol; j++)
            {
                if(max < FABS(MATRIX(inv, i, j)))
                {
                    max = FABS(MATRIX(inv, i, j));
                    is[k] = i;
                    js[k] = j;
                }
            }
        }

        if(max < 0.0001)
        {   //!     
            ret = -1;
            goto end;
        }

        //  
        if(is[k] != k)
        {
            matrix_swap_row(inv, k, is[k]);
        }
        if(js[k] != k)
        {
            matrix_swap_col(inv, k, js[k]);
        }

        MATRIX(inv, k, k) = 1 / MATRIX(inv, k, k);

        for(j=0; jcol; j++)
        {
            if(j != k)
                MATRIX(inv, k, j) *= MATRIX(inv, k, k);
        }
        for(i=0; irow; i++)
        {
            if(i != k)
            {
                for(j=0; jcol; j++)
                {
                    if(j != k)
                        MATRIX(inv, i, j) -= MATRIX(inv, i, k) * MATRIX(inv, k, j);
                }
            }
        }
        for(i=0; irow; i++)
        {
            if(i != k)
                MATRIX(inv, i, k) *= -MATRIX(inv, k, k);
        }

    }

    //  
    //   row  is[k], column  js[k]
    //   ：row  js[k], column  is[k]
    for(k=inv->row-1; k>=0; k--)
    {
        if(js[k] != k)
        {
            matrix_swap_row(inv, k, js[k]);
        }
        if(is[k] != k)
        {
            matrix_swap_col(inv, k, is[k]);
        }
    }

end:
    return ret;;
}

GitHub:すべてのコード
転載先:https://www.cnblogs.com/electron/p/5813789.html