function [solution,sigma,residual]=linsvd(matrix,y,plotarg);

%% Программа вычисляет коэффициенты линейной модели МНК
%% (метода наименьших квадратов) вида
%
% y(i) = a(1)*x(1) + ... +a(M)*x(M), i=1,2,...,N > M,
%
% которая может быть записана в матричном виде A * X = Y,
% где вектор-столбец X - неизвестные [x(1) x(2) ... x(M)]',
% а Y - вектор-столбец "левых частей" [y(1) y(2) ... y(N)]',
% матрица A (M столбцов, N > M строк) - матрица известных
% коэффициентов a(i,j) при неизвестных x(j), т.е.
%     |a(1,1)   a(1,2)   ... a(1,M)  |
%     |a(2,1)   a(2,2)   ... a(2,M)  |
% A = |..............................|
%     |a(N-1,1) a(N-1,2) ... a(N-1,M)|
%     |a(N,1)   a(N,2)   ... a(N,M)  |
%
% K - число неизвестных, N > M - число условных уравнений
%
% Примечание: по умолчанию данные считаются равноточными,
% но для неравноточных наблюдений/измерений следует
% перед вызовом программы поделить все коэффициенты a(i,j)
% i-й строки и левые части y(i) на ошибки левых частей sigma_y(i)
% В такой модификации программа реализует chi^2-минимизацию
%
% Обращение к программе:
%
% [solution,sigma[,residual]]=linsvd(matrix,y[,plotarg]);
%
% (в скобках указаны необязательные входные (выходные) параметры
%% Вход:
% matrix        - матрица A, столбцы которой являются коэффициентами
%                 при неизвестных для системы условных уравнений
% y             - вектор-столбец "левой части" y
% [plotarg]     - если задан, выводится график зависимости от него
%                 остаточных уклонений (residual)(необязательный
%                 параметр)
%% Выход:
% solution(:,1) - SVD МНК-решение линейной системы
% solution(:,2) - ошибки линейных коэффициентов
% sigma         - среднеквадратичная невязка (в смысле chi^2,
%                 т.е. средневзвешенная)
% [residual]    - массив уклонений от решения (необязательный параметр)
%
%% Автор: А.С.Расторгуев, физфак МГУ, ГАИШ МГУ, 2000-2013
%%

[n,m]=size(matrix);

%% SVD-разложение исходной матрицы:
[U,S,V]=svd(matrix,0);
for i=1:m
% Сингулярные числа:
   w(i)=S(i,i);
   S(i,i)=1/S(i,i);
end
solution=zeros(m,2);
% Решение  SVD:
solution(:,1)=V*S*U'*y;
zsol=matrix*solution(:,1);
%% Вычислим сумму квадратов невязок chi^2:
chi2=0;
residual(:,1)=y(:)-zsol(:);
for i=1:n
   chi2=chi2+residual(i)^2;
end
chi2=sqrt(chi2/(n-m));
sigma=chi2;
%% Ошибки неизвестных:
for j=1:m
   solution(j,2)=0;
   for i=1:m
      solution(j,2)=solution(j,2)+(V(j,i)*S(i,i))^2;
   end
   solution(j,2)=sqrt(solution(j,2))*chi2;
end
%
%% График остаточных уклонений
if nargin>2,
%    subplot(211)
%    plot(plotarg,zsol,'g.',plotarg,y,'r.');grid;xlabel('Plot argument');ylabel('Dependent var');
%    subplot(212)
    plot(plotarg,residual,'b.');grid;xlabel('Plot argument');ylabel('Residuals');
end;