Описание подпрограммы

Этот алгоритм является одним из немногих на сайте, которые требуют обратной коммуникации. Для его работы требуется операция умножения матрицы на вектор, которую подпрограмма IterativeEstimate1Norm не в состоянии осуществить самостоятельно - матрица A задана в неявной форме, и подпрограмма просто не может знать, как осуществлять умножение на неё.

Поэтому работа организована следующим образом: при первом запуске в подпрограмму передается параметр KASE, равный 0, что обозначает начало работы итерационного алгоритма. Если после возврата из подпрограммы параметр KASE неравен 0, то вызвавшая подпрограмма должна осуществить умножение возвращенного вектора X на матрицу A (если параметр равен 1) или на матрицу A ^T (если параметр равен 2) и повторно вызвать подпрограмму, оставив прочие параметры без изменений. Если после возврата из подпрограммы параметр KASE равен 0, то итерационный алгоритм завершил работу и подпрограмма вернула оценку нормы матрицы.

Более подпробно список параметров подпрограммы IterativeEstimate1Norm и работа с ней рассмотрены в комментариях. Имеет смысл обратить внимание на подпрограмму DemoIterativeEstimate1Norm, которая демонстрирует работу с подпрограммой IterativeEstimate1Norm.

(***********************************************************************Этот код сгенерирован транслятором AlgoPascal в рамках проекта ALGLIB Условия, на которых распространяется этот код, можно получить по адресу http://alglib.sources.ru/copyrules.php***********************************************************************) uses Math, Ap; procedure IterativeEstimate1Norm(N: Integer; var V: TReal1DArray; var X: TReal1DArray; var ISGN: TInteger1DArray; var EST: Double; var KASE: Integer); forward; function DemoIterativeEstimate1Norm(const A: TReal2DArray; N: Integer):Double; forward; (*************************************************************************Оценка нормы матрицы Алгоритм позволяет оценить 1-норму квадратной матрицы A при условии, чтодоступна операция умножения матрицы A на вектор (используется итеративныйметод). Алгоритм имеет смысл применять, если нет возможности сформироватьматрицу A (например, при оценке нормы обратной матрицы). Алгоритм использует обратную коммуникацию для умножения вектора наматрицу. Если при возврате из процедуры параметр KASE равен 0, процедурауспешно завершила работу, в противном случае процедура вернула управлениедля того, чтобы вызвавшая её программа осуществила операцию умножениявозвращенного вектора на матрицу A, и повторно вызвала подпрограмму. Подпрограмма DemoIterativeEstimateNorm демонстрирует это на простом примере. Параметры: N - размер задачи V - рабочий вектор. Самостоятельно инициализируется подпрограммой при первом вызове, передается в неё при повторных вызовах. X - при промежуточном возврате (KASE<>0) содержит вектор, который вызывающая подпрограмма должна заместить вектором A * X, если KASE=1, A^T * X, если KASE=2, и повторно вызвать подпрограмму, оставив остальные параметры без изменения. Массив с нумерацией элементов [1..N] ISGN- рабочий вектор. Самостоятельно инициализируется подпрограммой при первом вызове, передается в неё при повторных вызовах. EST - при финальном возврате (KASE=0) содержит нижнюю оценку нормы матрицы A. KASE- при первоначальном вызове подпрограммы параметр должен быть равен 0. При финальном возврате параметр равен 0 (EST содержит норму матрицы), при промежуточном возврате параметр равен 1 или 2, сообщая о том, какая операция должна быть проведена с вектором X. -- LAPACK auxiliary routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University February 29, 1992*************************************************************************) procedure IterativeEstimate1Norm(N: Integer; var V: TReal1DArray; var X: TReal1DArray; var ISGN: TInteger1DArray; var EST: Double; var KASE: Integer); var ITMAX: Integer; I: Integer; T: Double; Flg: Boolean; PosITER: Integer; PosJ: Integer; PosJLAST: Integer; PosJUMP: Integer; PosALTSGN: Integer; PosESTOLD: Integer; PosTEMP: Integer; i_: Integer; begin ITMAX:= 5; PosALTSGN:= N+1; PosESTOLD:= N+2; PosTEMP:= N+3; PosITER:= N+1; PosJ:= N+2; PosJLAST:= N+3; PosJUMP:= N+4; if KASE=0 then begin SetLength(V, N+3+1); SetLength(X, N+1); SetLength(ISGN, N+4+1); T:= 1/N; I:=1; while I<=N do begin X[I]:= T; Inc(I); end; KASE:= 1; ISGN[PosJUMP]:= 1; Exit; end; if ISGN[PosJUMP]=1 then begin if N=1 then begin V[1]:= X[1]; EST:= ABSReal(V[1]); KASE:= 0; Exit; end; EST:= 0; I:=1; while I<=N do begin EST:= EST+AbsReal(X[I]); Inc(I); end; I:=1; while I<=N do begin if X[I]>=0 then begin X[I]:= 1; end else begin X[I]:= -1; end; ISGN[I]:= Sign(X[I]); Inc(I); end; KASE:= 2; ISGN[PosJUMP]:= 2; Exit; end; if ISGN[PosJUMP]=2 then begin ISGN[PosJ]:= 1; I:=2; while I<=N do begin if AbsReal(X[I])>AbsReal(X[ISGN[PosJ]]) then begin ISGN[PosJ]:= I; end; Inc(I); end; ISGN[PosITER]:= 2; I:=1; while I<=N do begin X[I]:= 0; Inc(I); end; X[ISGN[PosJ]]:= 1; KASE:= 1; ISGN[PosJUMP]:= 3; Exit; end; if ISGN[PosJUMP]=3 then begin for i_:= 1 to N do begin V[i_]:= X[i_]; end; V[PosESTOLD]:= EST; EST:= 0; I:=1; while I<=N do begin EST:= EST+AbsReal(V[I]); Inc(I); end; Flg:= False; I:=1; while I<=N do begin if (X[I]>=0) and (ISGN[I]<0) or (X[I]<0) and (ISGN[I]>=0) then begin Flg:= True; end; Inc(I); end; if not Flg or (EST<=V[PosESTOLD]) then begin V[PosALTSGN]:= 1; I:=1; while I<=N do begin X[I]:= V[PosALTSGN]*(1+(I-1)/(N-1)); V[PosALTSGN]:= -V[PosALTSGN]; Inc(I); end; KASE:= 1; ISGN[PosJUMP]:= 5; Exit; end; I:=1; while I<=N do begin if X[I]>=0 then begin X[I]:= 1; ISGN[I]:= 1; end else begin X[I]:= -1; ISGN[I]:= -1; end; Inc(I); end; KASE:= 2; ISGN[PosJUMP]:= 4; Exit; end; if ISGN[PosJUMP]=4 then begin ISGN[PosJLAST]:= ISGN[PosJ]; ISGN[PosJ]:= 1; I:=2; while I<=N do begin if AbsReal(X[I])>AbsReal(X[ISGN[PosJ]]) then begin ISGN[PosJ]:= I; end; Inc(I); end; if (X[ISGN[PosJLAST]]<>ABSReal(X[ISGN[PosJ]])) and (ISGN[PosITER]<ITMAX) then begin ISGN[PosITER]:= ISGN[PosITER]+1; I:=1; while I<=N do begin X[I]:= 0; Inc(I); end; X[ISGN[PosJ]]:= 1; KASE:= 1; ISGN[PosJUMP]:= 3; Exit; end; V[PosALTSGN]:= 1; I:=1; while I<=N do begin X[I]:= V[PosALTSGN]*(1+(I-1)/(N-1)); V[PosALTSGN]:= -V[PosALTSGN]; Inc(I); end; KASE:= 1; ISGN[PosJUMP]:= 5; Exit; end; if ISGN[PosJUMP]=5 then begin V[PosTEMP]:= 0; I:=1; while I<=N do begin V[PosTEMP]:= V[PosTEMP]+AbsReal(X[I]); Inc(I); end; V[PosTEMP]:= 2*V[PosTEMP]/(3*N); if V[PosTEMP]>EST then begin for i_:= 1 to N do begin V[i_]:= X[i_]; end; EST:= V[PosTEMP]; end; KASE:= 0; Exit; end; end; (*************************************************************************Демонстрация работы с подпрограммой IterativeEstimateNorm. Подпрограмма принимает матрицу размером NxN (массив с нумерацией элементов[1..N, 1..N]) и возвращает норму, оцененную при помощи подпрограммы,демонстрируя принцип работы с подпрограммой IterativeEstimateNorm. -- ALGLIB -- Copyright 2005 by Bochkanov Sergey*************************************************************************) function DemoIterativeEstimate1Norm(const A: TReal2DArray; N: Integer):Double; var I: Integer; S: Double; X: TReal1DArray; T: TReal1DArray; V: TReal1DArray; IV: TInteger1DArray; KASE: Integer; i_: Integer; begin KASE:= 0; SetLength(T, N+1); IterativeEstimate1Norm(N, V, X, IV, Result, KASE); while KASE<>0 do begin if KASE=1 then begin I:=1; while I<=N do begin S:= 0; for i_:= 1 to N do begin S:= S + A[I,i_]*X[i_]; end; T[I]:= S; Inc(I); end; end else begin I:=1; while I<=N do begin S:= 0; for i_:= 1 to N do begin S:= S + A[i_,I]*X[i_]; end; T[I]:= S; Inc(I); end; end; for i_:= 1 to N do begin X[i_]:= T[i_]; end; IterativeEstimate1Norm(N, V, X, IV, Result, KASE); end; end;