Текст подпрограммы и версий
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) в произведение А = UDVT, где U, V - ортогональные матрицы порядка М и N соответственно, D - прямоугольная матрица размера M на N, у которой все элементы, кpоме диагональных элементов dk = dkk, k = 1, 2, ..., N, равны нулю.

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

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) = dk вычисленных сингулярных чисел матрицы 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 |