Приведение к Хессенбергу с помощью Гивенса

[Kirillko93]

Помогите разобраться! Мне нужно написать программу на $C++$ для приведения матрицы 6 на 6 к матрице Хессенберга с помощью преобразований Гивенса.
Не могу до конца понять, каким образом нужно в матрице $T$ расставлять значения $c$ и $s$. И что делать в последствии с матрицей $G$. Нужно ли ее куда-то записывать после очередного получения путем умножения $T_{kj}^T \cdot A \cdot T_{kj}$ и, если да, то куда именно? Сколько промежуточных матриц $G$ (т.е. сколько шагов) мы должны получить?Большое спасибо!

[AndrewN]

Kirillko93 в сообщении #589845 писал(а):

Помогите разобраться! Мне нужно написать программу на $C++$ для приведения матрицы 6 на 6 к матрице Хессенберга с помощью преобразований Гивенса.
Не могу до конца понять, каким образом нужно в матрице $T$ расставлять значения $c$ и $s$. И что делать в последствии с матрицей $G$. Нужно ли ее куда-то записывать после очередного получения путем умножения $T_{kj}^T \cdot A \cdot T_{kj}$ и, если да, то куда именно? Сколько промежуточных матриц $G$ (т.е. сколько шагов) мы должны получить?Большое спасибо!

Вращение обнуляет один элемент матрицы, а в форме Хессенберга нужно обнулить $N_0=(n-1)(n-2)/2$ элементов. Отсюда понятно число вращений.

Для обнуления элемента (индексы нумеруются, как массивы в С) $a_{20}$ составим матрицу:

$$R_{12} = \begin{pmatrix} c_{11} & -s_{12} \\ s_{21} & c_{22} \end{pmatrix} $$

которая в остальных местах единичная, а индексы обозначают расстановку коэффициентов $(c, s)$.

Коэффициенты $(c, s)$ вычисляются как $c=a_{10}/\sqrt{a_{10}^2+a_{20}^2}$ и $s=-a_{20}/\sqrt{a_{10}^2+a_{20}^2}$.

Остальные вращения вычисляются аналогично. Для каждого вращения берутся элементы активной подматрицы.

Хранить нужно $N_0$ пар коэффициентов $(c, s)$, т.е. $2N_0$ вещественных переменных.

А вообще, перед тем, как начинать что-то программировать, про это что-то хорошо бы и в учебнике почитать. Например, Х.Д.Икрамов, Численное решение матричных уравнений. Ортогональные методы. М., Наука, 1984, стр. 26-28. Обратите внимание на окончание - на стр. 28.

// ---------------------------------------------------------------------------
// Copyright (C) 2012 <andrew dot nesterov at techemail dot com>
// This code is not in public domain. Use it for educational purposes
// only.
// ---------------------------------------------------------------------------

// transform matrix to upper hessenberg form using givens rotations

#include <math.h>

#define N       5
#define sqr(x)  ((x)*(x))

// construct givens rotations R(k,i) to zero a[i][k-1]
// and possibly a[k-1][i], if matrix is symmetric
void crot (float akl, float ail, float *c, float *s)
{
    float b;

    b = 1.0e0f / sqrt (sqr (akl) + sqr (ail));
    *c =  akl * b;
    *s = -ail * b;

    return;
}

// apply left and right givens rotations R(k,i) to zero a[i][k-1]
// and possibly a[k-1][i], if matrix is symmetric
void grot (int n, int k, int i, float c, float s, float a[N][N])
{
    int j;
    float ai[N], ak[N];

    // left multiply (two dot products with two terms)
    for (j = 0; j < n; j++)                 // all cols
    {
        // k-th row
        ak[j] = c * a[k][j] - s * a[i][j];
        // i-th row
        ai[j] = s * a[k][j] + c * a[i][j];
    }

    // copy temp arrays
    for (j = 0; j < n; j++)                 // all cols
    {
        // k-th row
        a[k][j] = ak[j];
        // i-th row
        a[i][j] = ai[j];
    }

    // right multiply (two dot products with two terms)
    for (j = 0; j < n; j++)                 // all rows
    {
        // k-th col
        ak[j] = c * a[j][k] - s * a[j][i];
        // i-th col
        ai[j] = s * a[j][k] + c * a[j][i];
    }

    // copy temp arrays
    for (j = 0; j < n; j++)                 // all rows
    {
        // k-th col
        a[j][k] = ak[j];
        // i-th col
        a[j][i] = ai[j];
    }

    return;
}

#include <stdio.h>

// debug output
void printm (int n, float a[N][N])
{
    int i, k;

    for (i = 0; i < n; i++)                 // row cycle
    {
        for (k = 0; k < n; k++)             // col cycle
        {
            printf ("%14.7f ", a[i][k]);
        }

        printf ("\n");
    }

    printf ("\n");

    return;
}

// transform matrix to upper hessenberg form
void hessen (int n, float a[N][N], float *cs)
{
    int i, j, k, l;
    float c, s, b;

    j = 0;

    for (k = 1; k < n - 1; k++)             // col cycle
    {
        for (i = k + 1; i < n; i++)         // row cycle
        {
            // construct givens rotations R(k,i) to zero a[i][k-1]
            // and a[k-1][i] if matrix is symmetric
            l = k - 1;
            crot (a[k][l], a[i][l], &c, &s);

            // store givens rotation (i,k)
            cs[2*j+0] = c;
            cs[2*j+1] = s;
            j += 1;

            // apply givens rotations R(k,i) to zero a[i][k-1]
            // if A is symmetric, a[k-1][i] is zeroed as well
            grot (n, k, i, c, s, a);

            // printm (n, a);                  // debug phase
        }
    }

    return;
}

// test givens rotations

#include <stdlib.h>

// initialize generic matrix
void minit (int n, float a[N][N])
{
    int i, k;

    // init generic matrix
    for (i = 0; i < n; i++)                 // row cycle
    {
        for (k = 0; k < n; k++)             // col cycle
        {
            a[i][k] = (float) rand () / (float) RAND_MAX;
        }
    }

    return;
}

// initialize symmetric matrix
void sinit (int n, float a[N][N])
{
    int i, k;

    // init symmetric matrix
    for (i = 0; i < n; i++)                 // row cycle
    {
        for (k = i; k < n; k++)             // col cycle
        {
            a[i][k] = (float) rand () / (float) RAND_MAX;
            if (i < k) a[k][i] = a[i][k];
        }
    }

    return;
}

// global data
float A[N][N];
// storage for up to 10 x 10 matrix transform
float cs[9*8/2];

void main (void)
{
    int i, k, n;

    n = N;

    minit  (n, A);                          // init generic matrix
    // sinit  (n, A);                          // init symmetric matrix

    printm (n, A);

    hessen (n, A, cs);

    printm (n, A);

    exit (0);
}