БЧА НИВЦ МГУ. ACP1R_C. Сингулярное разложение матриц общего вида

Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.su/lib_na/cat/ac_htm_c/acp1r_c.htm
Дата изменения: Tue Apr 28 11:43:14 2015
Дата индексирования: Sun Apr 10 02:05:32 2016
Кодировка: Windows-1251

Текст подпрограммы и версий
acp1r_c.zip , acp1d_c.zip

Тексты тестовых примеров
tacp1r_c.zip , tacp1d_c.zip

Подпрограмма: acp1r_c

Назначение

Сингулярное разложение прямоугольной матрицы размера m на n (m ≥ n).

Математическое описание

Подпрограмма acp1r_c осуществляет разложение прямоугольной матрицы A размера M на N (M ≥ N) в произведение А = UDV^T, где U, V - ортогональные матрицы порядка М и N соответственно, D - прямоугольная матрица размера M на N, у которой все элементы, кpоме диагональных элементов d_k = d_kk, k = 1, 2, ..., N, равны нулю.

Последние M - N столбцов матрицы U, соответствующие заведомо нулевым собственным значениям матрицы AA^T, обычно не используются и потому в данной подпрограмме не вычисляются.

Forsythe G.E., Malcolm M.A., Moler C.B. Computer Methods for Mathematical Computation. Prentice-Hall, 1977.

Использование

    int acp1r_c (integer *nm, integer *m, integer *n, real *a,
                 real *w, logical *matu, real *u, logical *matv, real *v, real *rv1, 
                 integer *ierr)

   Параметры

nm -	число стpок в двумерных массивах a, u, v, причем nm ≥ m (тип: целый);
m -	число стpок матриц a и u (тип: целый);
n -	число столбцов матриц a и u и порядок матрицы v (тип: целый);
a -	вещественный двумерный массив размера nm на n, в первых m стpоках которого задается исходная матрица;
w -	вещественный вектоp длины n значений w(k) = d_k вычисленных сингулярных чисел матрицы a;
matu -	признак необходимости вычисления матрицы u: при matu = TRUE_ матрица u вычисляется, при matu = FALSE_ - нет (тип: логический);
u -	вещественный двумерный массив размера nm на n, в котоpом при matu = TRUE_ содержатся вычисленные первые n столбцов матрицы u; при matu = FALSE_ массив u используется как рабочий;
matv -	признак необходимости вычисления матрицы v: при matv = TRUE_ матрица v вычисляется, при matv = FALSE_ - нет (тип: логический);
v -	вещественный двумерный массив размера nm на n, в котоpом при matv = TRUE_ содержится вычисленная матрица v; при matv = FALSE_ массив v не используется;
r -	вещественный рабочий вектоp длины n;
ierr -	целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:

ierr = k, если для вычисления k - го сингулярного числа потребовалось более 30 итераций; в этом случае правильно вычислены лишь сингулярные числа и сингулярные векторы с индексами k+1, k+2, ..., n.

   Версии

acp1d_c -

сингулярное разложение прямоугольной матрицы размера m на n (m ≥ n), заданной с удвоенной точностью.

   Вызываемые подпрограммы

utac10_c -

подпрограмма выдачи диагностических сообщений при работе подпрограмм acp1r_c и acp1d_c.

   Замечания по использованию

	1.	При обращении к подпрограмме acp1d_c параметры a, w, u, v, r должны иметь тип double.
	2.	B результате работы подпрограммы массив a сохраняется.
	3.	Допустимо совпадение параметров u или v с параметром a.

Пример использования

int main(void)
{
    /* System generated locals */
    int i__1, i__2;

    /* Local variables */
    static int ierr;
    static float work[5];
    extern int acp1r_c(int *, int *, int *, float *, float *,
                      logical *, float *, logical *, float *, float *, int *);
    static float a[15]   /* was [5][3] */;
    static int i__, j, m, n;
    static float u[15] /* was [5][3] */,
                 v[15] /* was [5][3] */, sigma[3];
    static int nm;

#define a_ref(a_1,a_2) a[(a_2)*5 + a_1 - 6]
#define u_ref(a_1,a_2) u[(a_2)*5 + a_1 - 6]
#define v_ref(a_1,a_2) v[(a_2)*5 + a_1 - 6]

    nm = 5;
    m = 5;
    n = 3;
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            a_ref(i__, j) = (float) (i__ + (j - 1) * m);
/* l1: */
        }
    }
    acp1r_c(&nm, &m, &n, a, sigma, &c_true, u, &c_true, v, work, &ierr);

    if (ierr != 0) {
        printf("\n trouble, ierr =  %5i \n",ierr);
    }
        printf("\n %16.7e %16.7e %16.7e \n\n",
            sigma[0], sigma[1], sigma[2]);
/* l3: */
    for (j = 1; j <= 5; ++j) {
        printf("\n  %16.7e %16.7e %16.7e \n",
            u_ref(j, 1), u_ref(j, 2), u_ref(j, 3));
/* l4: */
    }
    for (j = 1; j <= 3; ++j) {
        printf("\n  %16.7e %16.7e %16.7e \n",
            v_ref(j, 1), v_ref(j, 2), v_ref(j, 3));
/* l5: */
    }
    return 0;
} /* main */


Результаты:
          ierr  =  0, 

               | 35.127         |
 sigma  =  |   2.465         |     , 
               |   1.587e-11 |

              | -0.202    0.890    0.408  | 
 v_ref  =  | -0.517    0.257   -0.816 |     , 
              | -0.832   -0.376    0.408 |

              | -0.355   -0.689    0.602 |
              | -0.399   -0.376   -0.517 |
 u_ref  =  | -0.443   -0.062   -0.306 |
              | -0.487    0.251   -0.244 |
              | -0.531    0.564    0.466 |