Mercurial > hg > octave-terminal
changeset 7132:b01db194c526
[project @ 2007-11-08 16:17:34 by jwe]
author | jwe |
---|---|
date | Thu, 08 Nov 2007 16:17:34 +0000 |
parents | a184bc985c37 |
children | 1d0d7be2d0f8 |
files | scripts/ChangeLog scripts/control/util/__outlist__.m scripts/control/util/__zgpbal__.m scripts/control/util/axis2dlim.m scripts/control/util/prompt.m scripts/control/util/sortcom.m scripts/control/util/zgfmul.m scripts/control/util/zgfslv.m scripts/control/util/zginit.m scripts/control/util/zgreduce.m scripts/control/util/zgrownorm.m scripts/control/util/zgscal.m |
diffstat | 12 files changed, 222 insertions(+), 177 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -3,7 +3,12 @@ * control/hinf/h2syn.m, control/hinf/hinf_ctr.m, control/hinf/hinfnorm.m, control/hinf/hinfsyn.m, control/hinf/hinfsyn_chk.m, control/hinf/is_dgkf.m, - control/hinf/wgt1o.m: Style fixes. + control/hinf/wgt1o.m, control/util/__outlist__.m, + control/util/__zgpbal__.m, control/util/axis2dlim.m, + control/util/prompt.m, control/util/sortcom.m, + control/util/zgfmul.m, control/util/zgfslv.m, + control/util/zginit.m, control/util/zgreduce.m, + control/util/zgrownorm.m, control/util/zgscal.m: Style fixes. 2007-11-08 David Bateman <dbateman@free.fr>
--- a/scripts/control/util/__outlist__.m +++ b/scripts/control/util/__outlist__.m @@ -49,33 +49,33 @@ function str_val = __outlist__ (name_list, tabchar, yd, ilist) - if( nargin < 1 | nargin > 4 ) + if (nargin < 1 || nargin > 4) print_usage (); endif - m = length(name_list); - if(nargin < 4) + m = length (name_list); + if (nargin < 4) ilist = 1:m; endif - if(nargin ==1) + if (nargin == 1) tabchar = ""; endif - if(nargin < 3) - yd = zeros(1,m); - elseif(isempty(yd)) - yd = zeros(1,m); + if (nargin < 3) + yd = zeros (1, m); + elseif (isempty (yd)) + yd = zeros (1, m); endif str_val = ""; - dstr = {""," (discrete)"}; - if((m >= 1) && (iscell(name_list))) - for ii=1:m - str_val = sprintf("%s%s%d: %s%s\n",str_val,tabchar, ilist(ii), ... - name_list{ii},dstr{yd(ii)+1}); + dstr = {"", " (discrete)"}; + if (m >= 1 && iscell (name_list)) + for ii = 1:m + str_val = sprintf ("%s%s%d: %s%s\n", str_val, tabchar, ilist(ii), + name_list{ii}, dstr{yd(ii)+1}); endfor else - str_val = sprintf("%sNone",tabchar); + str_val = sprintf ("%sNone", tabchar); endif endfunction
--- a/scripts/control/util/__zgpbal__.m +++ b/scripts/control/util/__zgpbal__.m @@ -47,67 +47,67 @@ function retsys = __zgpbal__ (Asys) - if( (nargin != 1) | (!isstruct(Asys))) + if (nargin != 1 || ! isstruct (Asys)) print_usage (); endif - Asys = sysupdate(Asys,"ss"); - [a,b,c,d] = sys2ss(Asys); + Asys = sysupdate (Asys, "ss"); + [a, b, c, d] = sys2ss (Asys); - [nn,mm,pp] = abcddim(a,b,c,d); + [nn, mm, pp] = abcddim (a, b, c, d); np1 = nn+1; nmp = nn+mm+pp; ## set up log vector zz, incidence matrix ff - zz = zginit(a,b,c,d); + zz = zginit (a, b, c, d); ## disp("__zgpbal__: zginit returns") ## zz ## disp("/__zgpbal__") - if (norm(zz)) + if (norm (zz)) ## generalized conjugate gradient approach - xx = zgscal(a,b,c,d,zz,nn,mm,pp); + xx = zgscal (a, b, c, d, zz, nn, mm, pp); - for i=1:nmp - xx(i) = floor(xx(i)+0.5); + for i = 1:nmp + xx(i) = floor (xx(i)+0.5); xx(i) = 2.0^xx(i); endfor ## now scale a ## block 1: a = sigma a inv(sigma) - for i=1:nn + for i = 1:nn a(i,1:nn) = a(i,1:nn)*xx(i); a(1:nn,i) = a(1:nn,i)/xx(i); endfor ## block 2: b= sigma a phi - for j=1:mm + for j = 1:mm j1 = j+nn; b(1:nn,j) = b(1:nn,j)*xx(j1); endfor - for i=1:nn + for i = 1:nn b(i,1:mm) = b(i,1:mm)*xx(i); endfor - for i=1:pp + for i = 1:pp i1 = i+nn+mm; ## block 3: c = psi C inv(sigma) c(i,1:nn) = c(i,1:nn)*xx(i1); endfor - for j=1:nn + for j = 1:nn c(1:pp,j) = c(1:pp,j)/xx(j); endfor ## block 4: d = psi D phi - for j=1:mm + for j = 1:mm j1 = j+nn; d(1:pp,j) = d(1:pp,j)*xx(j1); endfor - for i=1:pp + for i = 1:pp i1 = i + nn + mm; d(i,1:mm) = d(i,1:mm)*xx(i1); endfor endif - retsys = ss(a,b,c,d); + retsys = ss (a, b, c, d); + endfunction -
--- a/scripts/control/util/axis2dlim.m +++ b/scripts/control/util/axis2dlim.m @@ -49,24 +49,24 @@ delv = (maxv-minv)/2; # breadth of the plot midv = (minv + maxv)/2; # midpoint of the plot axmid = [midv(1), midv(1), midv(2), midv(2)]; - axdel = [-0.1, 0.1,-0.1,0.1]; # default plot width (if less than 2-d data) - if(max(delv) == 0) - if(midv(1) != 0) - axdel(1:2) = [-0.1*midv(1),0.1*midv(1)]; + axdel = [-0.1, 0.1, -0.1, 0.1]; # default plot width (if less than 2-d data) + if (max (delv) == 0) + if (midv(1) != 0) + axdel(1:2) = [-0.1*midv(1), 0.1*midv(1)]; endif - if(midv(2) != 0) - axdel(3:4) = [-0.1*midv(2),0.1*midv(2)]; + if (midv(2) != 0) + axdel(3:4) = [-0.1*midv(2), 0.1*midv(2)]; endif else ## they're at least one-dimensional tolv = max(1e-8, 1e-8*abs(midv)); - if(abs(delv(1)) >= tolv(1)) + if (abs (delv(1)) >= tolv(1)) axdel(1:2) = 1.1*[-delv(1),delv(1)]; endif - if(abs(delv(2)) >= tolv(2)) + if (abs (delv(2)) >= tolv(2)) axdel(3:4) = 1.1*[-delv(2),delv(2)]; endif endif axvec = axmid + axdel; + endfunction -
--- a/scripts/control/util/prompt.m +++ b/scripts/control/util/prompt.m @@ -41,7 +41,7 @@ print_usage (); elseif (nargin == 0) str = "\n ---- Press a key to continue ---"; - elseif (! ischar (str) ) + elseif (! ischar (str)) error ("prompt: input must be a string"); endif
--- a/scripts/control/util/sortcom.m +++ b/scripts/control/util/sortcom.m @@ -53,42 +53,47 @@ function [yy, idx] = sortcom (xx, opt) - if( nargin < 1 | nargin > 2 ) + if (nargin < 1 || nargin > 2) print_usage (); - elseif( !(isvector(xx) | isempty(xx) )) - error("sortcom: first argument must be a vector"); + elseif (! (isvector (xx) || isempty (xx))) + error ("sortcom: first argument must be a vector"); endif - if(nargin == 1) opt = "re"; + if (nargin == 1) + opt = "re"; else - if (!ischar(opt)) - error("sortcom: second argument must be a string"); + if (! ischar (opt)) + error ("sortcom: second argument must be a string"); endif endif - if(isempty(xx)) + if (isempty (xx)) yy = idx = []; else - if(strcmp(opt,"re")) datavec = real(xx); - elseif(strcmp(opt,"im")) datavec = imag(xx); - elseif(strcmp(opt,"mag")) datavec = abs(xx); - else error(["sortcom: invalid option = ", opt]) + if (strcmp (opt, "re")) + datavec = real (xx); + elseif (strcmp (opt, "im")) + datavec = imag (xx); + elseif (strcmp (opt, "mag")) + datavec = abs (xx); + else + error ("sortcom: invalid option = %s", opt); endif - [datavec,idx] = sort(datavec); + [datavec, idx] = sort (datavec); yy= xx(idx); - if(strcmp(opt,"re") | strcmp(opt,"mag")) + if (strcmp (opt, "re") || strcmp (opt, "mag")) ## sort so that complex conjugate pairs appear together - ddiff = diff(datavec); - zidx = find(ddiff == 0); + ddiff = diff (datavec); + zidx = find (ddiff == 0); ## sort common datavec values - if(!isempty(zidx)) - for iv=create_set(datavec(zidx)) - vidx = find(datavec == iv); - [vals,imidx] = sort(imag(yy(vidx))); + if (! isempty (zidx)) + for iv = create_set (datavec(zidx)) + vidx = find (datavec == iv); + [vals, imidx] = sort (imag (yy(vidx))); yy(vidx) = yy(vidx(imidx)); idx(vidx) = idx(vidx(imidx)); endfor
--- a/scripts/control/util/zgfmul.m +++ b/scripts/control/util/zgfmul.m @@ -37,48 +37,67 @@ print_usage (); endif - [n,m] = size(b); - [p,m1] = size(c); + [n, m] = size (b); + [p, m1] = size (c); nm = n+m; - y = zeros(nm+p,1); + y = zeros (nm+p, 1); ## construct F column by column - for jj=1:n - Fj = zeros(nm+p,1); + for jj = 1:n + Fj = zeros (nm+p, 1); ## rows 1:n: F1 - aridx = complement(jj,find(a(jj,:) != 0)); - acidx = complement(jj,find(a(:,jj) != 0)); - bidx = find(b(jj,:) != 0); - cidx = find(c(:,jj) != 0); + aridx = complement (jj, find (a(jj,:) != 0)); + acidx = complement (jj, find (a(:,jj) != 0)); + bidx = find (b(jj,:) != 0); + cidx = find (c(:,jj) != 0); Fj(aridx) = Fj(aridx) - 1; # off diagonal entries of F1 Fj(acidx) = Fj(acidx) - 1; ## diagonal entry of F1 - Fj(jj) = length(aridx)+length(acidx) + length(bidx) + length(cidx); + Fj(jj) = length (aridx) + length (acidx) + length (bidx) + length (cidx); - if(!isempty(bidx)) Fj(n+bidx) = 1; endif # B' incidence - if(!isempty(cidx)) Fj(n+m+cidx) = -1; endif # -C incidence + ## B' incidence + if (! isempty (bidx)) + Fj(n+bidx) = 1; + endif + + ## -C incidence + if (! isempty (cidx)) + Fj(n+m+cidx) = -1; + endif y = y + x(jj)*Fj; # multiply by corresponding entry of x endfor - for jj=1:m - Fj = zeros(nm+p,1); - bidx = find(b(:,jj) != 0); - if(!isempty(bidx)) Fj(bidx) = 1; endif # B incidence - didx = find(d(:,jj) != 0); - if(!isempty(didx)) Fj(n+m+didx) = 1; endif # D incidence + for jj = 1:m + Fj = zeros (nm+p, 1); + bidx = find (b(:,jj) != 0); + ## B incidence + if (! isempty (bidx)) + Fj(bidx) = 1; + endif + didx = find (d(:,jj) != 0); + ## D incidence + if (! isempty (didx)) + Fj(n+m+didx) = 1; + endif Fj(n+jj) = length(bidx) + length(didx); # F2 is diagonal y = y + x(n+jj)*Fj; # multiply by corresponding entry of x endfor - for jj=1:p - Fj = zeros(nm+p,1); - cidx = find(c(jj,:) != 0); - if(!isempty(cidx)) Fj(cidx) = -1; endif # -C' incidence + for jj = 1:p + Fj = zeros (nm+p, 1); + cidx = find (c(jj,:) != 0); + ## -C' incidence + if (! isempty (cidx)) + Fj(cidx) = -1; + endif didx = find(d(jj,:) != 0); - if(!isempty(didx)) Fj(n+didx) = 1; endif # D' incidence - Fj(n+m+jj) = length(cidx) + length(didx); # F2 is diagonal + ## D' incidence + if (! isempty (didx)) + Fj(n+didx) = 1; + endif + Fj(n+m+jj) = length (cidx) + length (didx); # F2 is diagonal y = y + x(n+m+jj)*Fj; # multiply by corresponding entry of x endfor
--- a/scripts/control/util/zgfslv.m +++ b/scripts/control/util/zgfslv.m @@ -32,41 +32,51 @@ endif nmp = n+m+p; - gam1 = (2*n)+m+p; gam2 = n+p; gam3 = n+m; + gam1 = (2*n)+m+p; + gam2 = n+p; + gam3 = n+m; - G1 = givens(sqrt(m),-sqrt(p))'; - G2 = givens(m+p,sqrt(n*(m+p)))'; + G1 = givens (sqrt (m), -sqrt (p))'; + G2 = givens (m+p, sqrt (n*(m+p)))'; x = b; ## 1) U1 e^n = sqrt(n)e_1^n ## 2) U2 e^m = sqrt(m)e_1^m ## 3) U3 e^p = sqrt(p)e_1^p - xdx1 = 1:n; xdx2 = n+(1:m); xdx3 = n+m+(1:p); - x(xdx1,1) = zgshsr(x(xdx1,1)); - x(xdx2,1) = zgshsr(x(xdx2,1)); - x(xdx3,1) = zgshsr(x(xdx3,1)); + xdx1 = 1:n; + xdx2 = n+(1:m); + xdx3 = n+m+(1:p); + + x(xdx1,1) = zgshsr (x(xdx1,1)); + x(xdx2,1) = zgshsr (x(xdx2,1)); + x(xdx3,1) = zgshsr (x(xdx3,1)); ## 4) Givens rotations to reduce stray non-zero elements - idx1 = [n+1,n+m+1]; idx2 = [1,n+1]; + idx1 = [n+1, n+m+1]; + idx2 = [1, n+1]; + x(idx1) = G1'*x(idx1); x(idx2) = G2'*x(idx2); ## 6) Scale x, then back-transform to get x - en = ones(n,1); em = ones(m,1); ep = ones(p,1); - lam = [gam1*en;gam2*em;gam3*ep]; + en = ones (n, 1); + em = ones (m, 1); + ep = ones (p, 1); + lam = [gam1*en; gam2*em; gam3*ep]; lam(1) = n+m+p; lam(n+1) = 1; # dummy value to avoid divide by zero - lam(n+m+1)=n+m+p; + lam(n+m+1) = n+m+p; - x = x ./ lam; x(n+1) = 0; # minimum norm solution + x = x ./ lam; + x(n+1) = 0; # minimum norm solution ## back transform now. x(idx2) = G2*x(idx2); x(idx1) = G1*x(idx1); - x(xdx3,1) = zgshsr(x(xdx3,1)); - x(xdx2,1) = zgshsr(x(xdx2,1)); - x(xdx1,1) = zgshsr(x(xdx1,1)); + x(xdx3,1) = zgshsr (x(xdx3,1)); + x(xdx2,1) = zgshsr (x(xdx2,1)); + x(xdx1,1) = zgshsr (x(xdx1,1)); endfunction
--- a/scripts/control/util/zginit.m +++ b/scripts/control/util/zginit.m @@ -38,58 +38,61 @@ print_usage (); endif - [nn,mm] = size(b); - [pp,mm] = size(d); + [nn, mm] = size (b); + [pp, mm] = size (d); nmp = nn+mm+pp; ## set up log vector zz - zz = zeros(nmp,1); + zz = zeros (nmp, 1); ## zz part 1: - for i=1:nn + for i = 1:nn ## nonzero off diagonal entries of a - if(nn > 1) - nidx = complement(i,1:nn); - a_row_i = a(i,nidx); a_col_i = a(nidx,i); - arnz = a_row_i(find(a_row_i != 0)); acnz = a_col_i(find(a_col_i != 0)); + if (nn > 1) + nidx = complement (i, 1:nn); + a_row_i = a(i,nidx); + a_col_i = a(nidx,i); + arnz = a_row_i(find (a_row_i != 0)); + acnz = a_col_i(find (a_col_i != 0)); else arnz = acnz = []; endif ## row of b - bidx = find(b(i,:) != 0); + bidx = find (b(i,:) != 0); b_row_i = b(i,bidx); ## column of c - cidx = find(c(:,i) != 0); + cidx = find (c(:,i) != 0); c_col_i = c(cidx,i); ## sum the entries - zz(i) = sum(log(abs(acnz))) - sum(log(abs(arnz))) ... - - sum(log(abs(b_row_i))) + sum(log(abs(c_col_i))); + zz(i) = sum (log (abs (acnz))) - sum (log (abs (arnz))) ... + - sum (log (abs (b_row_i))) + sum (log (abs (c_col_i))); endfor ## zz part 2: - bd = [b;d]; - for i=1:mm + bd = [b; d]; + for i = 1:mm i1 = i+nn; ## column of [b;d] - bdidx = find(bd(:,i) != 0); + bdidx = find (bd(:,i) != 0); bd_col_i = bd(bdidx,i); - zz(i1) = sum(log(abs(bd_col_i))); + zz(i1) = sum (log (abs(bd_col_i))); endfor ## zz part 3: cd = [c, d]; - for i=1:pp + for i = 1:pp i1 = i+nn+mm; - cdidx = find(cd(i,:) != 0); + cdidx = find (cd(i,:) != 0); cd_row_i = cd(i,cdidx); - zz(i1) = -sum(log(abs(cd_row_i))); + zz(i1) = -sum (log (abs (cd_row_i))); endfor ## now set zz as log base 2 - zz = zz*(1/log(2)); + zz *= (1 / log (2)); + endfunction
--- a/scripts/control/util/zgreduce.m +++ b/scripts/control/util/zgreduce.m @@ -31,43 +31,43 @@ ## SYS_INTERNAL accesses members of system data structure - is_digital(Asys); # make sure it's pure digital/continuous + is_digital (Asys); # make sure it's pure digital/continuous exit_1 = 0; # exit_1 = 1 or 2 on exit of loop - if(Asys.n + Asys.nz == 0) + if (Asys.n + Asys.nz == 0) exit_1 = 2; # there are no finite zeros endif while (! exit_1) - [Q,R,Pi] = qr(Asys.d); # compress rows of D + [Q, R, Pi] = qr (Asys.d); # compress rows of D Asys.d = Q'*Asys.d; Asys.c = Q'*Asys.c; ## check row norms of Asys.d - [sig,tau] = zgrownorm(Asys.d,meps); + [sig, tau] = zgrownorm (Asys.d, meps); ## disp("=======================================") ## disp(["zgreduce: meps=",num2str(meps), ", sig=",num2str(sig), ... ## ", tau=",num2str(tau)]) ## sysout(Asys) - if(tau == 0) + if (tau == 0) exit_1 = 1; # exit_1 - reduction complete and correct else Cb = Db = []; - if(sig) + if (sig) Cb = Asys.c(1:sig,:); Db = Asys.d(1:sig,:); endif - Ct =Asys.c(sig+(1:tau),:); + Ct = Asys.c(sig+(1:tau),:); ## compress columns of Ct - [pp,nn] = size(Ct); + [pp, nn] = size (Ct); rvec = nn:-1:1; - [V,Sj,Pi] = qr(Ct'); + [V, Sj, Pi] = qr (Ct'); V = V(:,rvec); - [rho,gnu] = zgrownorm(Sj,meps); + [rho, gnu] = zgrownorm (Sj, meps); ## disp(["zgreduce: rho=",num2str(rho),", gnu=",num2str(gnu)]) ## Cb @@ -75,22 +75,22 @@ ## Ct ## Sj' - if(rho == 0) + if (rho == 0) exit_1 = 1; # exit_1 - reduction complete and correct - elseif(gnu == 0) + elseif (gnu == 0) exit_1 = 2; # there are no zeros at all else mu = rho + sig; ## update system with Q - M = [Asys.a , Asys.b ]; - [nn,mm] = size(Asys.b); + M = [Asys.a, Asys.b ]; + [nn, mm] = size (Asys.b); pp = rows(Asys.d); - Vm =[V,zeros(nn,mm) ; zeros(mm,nn), eye(mm)]; - if(sig) + Vm =[V, zeros(nn,mm); zeros(mm,nn), eye(mm)]; + if (sig) M = [M; Cb, Db]; - Vs =[V',zeros(nn,sig) ; zeros(sig,nn), eye(sig)]; + Vs =[V', zeros(nn,sig); zeros(sig,nn), eye(sig)]; else Vs = V'; endif @@ -130,16 +130,16 @@ ## disp(["zgreduce: while loop done: exit_1=",num2str(exit_1)]); - if(exit_1 == 2) + if (exit_1 == 2) ## there are no zeros at all! Asys.a = Asys.b = Asys.c = []; endif ## update dimensions - if(is_digital(Asys)) - Asys.nz = rows(Asys.a); + if (is_digital (Asys)) + Asys.nz = rows (Asys.a); else - Asys.n = rows(Asys.a); + Asys.n = rows (Asys.a); endif retsys = Asys;
--- a/scripts/control/util/zgrownorm.m +++ b/scripts/control/util/zgrownorm.m @@ -38,4 +38,3 @@ tau = sum (rownorm <= meps); endfunction -
--- a/scripts/control/util/zgscal.m +++ b/scripts/control/util/zgscal.m @@ -44,27 +44,27 @@ nmp = n+m+p; ## x_0 = x_{-1} = 0, r_0 = z - x = zeros(nmp,1); + x = zeros (nmp, 1); xk1 = x; xk2 = x; rk1 = z; k = 0; ## construct balancing least squares problem - F = eye(nmp); - for kk=1:nmp - F(1:nmp,kk) = zgfmul(a,b,c,d,F(:,kk)); + F = eye (nmp); + for kk = 1:nmp + F(1:nmp,kk) = zgfmul (a, b, c, d, F(:,kk)); endfor - [U,H,k1] = krylov(F,z,nmp,1e-12,1); - if(!issquare(H)) - if(columns(H) != k1) - error("zgscal(tzero): k1=%d, columns(H)=%d",k1,columns(H)); - elseif(rows(H) != k1+1) - error("zgscal: k1=%d, rows(H) = %d",k1,rows(H)); - elseif ( norm(H(k1+1,:)) > 1e-12*norm(H,"inf") ) + [U, H, k1] = krylov (F, z, nmp, 1e-12, 1); + if (! issquare (H)) + if (columns (H) != k1) + error ("zgscal(tzero): k1=%d, columns(H)=%d", k1, columns (H)); + elseif (rows (H) != k1+1) + error ("zgscal: k1=%d, rows(H) = %d", k1, rows (H)); + elseif (norm (H(k1+1,:)) > 1e-12*norm (H, "inf")) zgscal_last_row_of_H = H(k1+1,:) - error("zgscal: last row of H nonzero (norm(H)=%e)",norm(H,"inf")) + error ("zgscal: last row of H nonzero (norm(H)=%e)", norm (H, "inf")) endif H = H(1:k1,1:k1); U = U(:,1:k1); @@ -72,8 +72,8 @@ ## tridiagonal H can still be rank deficient, so do permuted qr ## factorization - [qq,rr,pp] = qr(H); # H = qq*rr*pp' - nn = rank(rr); + [qq, rr, pp] = qr (H); # H = qq*rr*pp' + nn = rank (rr); qq = qq(:,1:nn); rr = rr(1:nn,:); # rr may not be square, but "\" does least xx = U*pp*(rr\qq'*(U'*z)); # squares solution, so this works @@ -86,38 +86,41 @@ ## so for now I'm solving it with the krylov routine. ## initialize residual error norm - rnorm = norm(rk1,1); + rnorm = norm (rk1, 1); xnorm = 0; - fnorm = 1e-12 * norm([a,b;c,d],1); + fnorm = 1e-12 * norm ([a, b; c, d], 1); - ## dummy defines for MATHTOOLS compiler - gamk2 = 0; omega1 = 0; ztmz2 = 0; + gamk2 = 0; + omega1 = 0; + ztmz2 = 0; ## do until small changes to x len_x = length(x); - while ((k < 2*len_x) & (xnorm> 0.5) & (rnorm>fnorm))|(k == 0) - k = k+1; + while ((k < 2*len_x && xnorm > 0.5 && rnorm > fnorm) || k == 0) + k++; ## solve F_d z_{k-1} = r_{k-1} - zk1= zgfslv(n,m,p,rk1); + zk1= zgfslv (n, m, p, rk1); ## Generalized CG iteration ## gamk1 = (zk1'*F_d*zk1)/(zk1'*F*zk1); ztMz1 = zk1'*rk1; - gamk1 = ztMz1/(zk1'*zgfmul(a,b,c,d,zk1)); + gamk1 = ztMz1/(zk1'*zgfmul (a, b, c, d, zk1)); - if(rem(k,len_x) == 1) omega = 1; - else omega = 1/(1-gamk1*ztMz1/(gamk2*omega1*ztmz2)); + if (rem (k, len_x) == 1) + omega = 1; + else + omega = 1/(1-gamk1*ztMz1/(gamk2*omega1*ztmz2)); endif ## store x in xk2 to save space xk2 = xk2 + omega*(gamk1*zk1 + xk1 - xk2); ## compute new residual error: rk = z - F xk, check end conditions - rk1 = z - zgfmul(a,b,c,d,xk2); - rnorm = norm(rk1); - xnorm = max(abs(xk1 - xk2)); + rk1 = z - zgfmul (a, b, c, d, xk2); + rnorm = norm (rk1); + xnorm = max (abs (xk1 - xk2)); ## printf("zgscal: k=%d, gamk1=%e, gamk2=%e, \nztMz1=%e ztmz2=%e\n", ... ## k,gamk1, gamk2, ztMz1, ztmz2); @@ -129,21 +132,22 @@ gamk2 = gamk1; omega1 = omega; ztmz2 = ztMz1; - [xk1,xk2] = swap(xk1,xk2); + [xk1, xk2] = swap (xk1, xk2); endwhile x = xk2; ## check convergence - if (xnorm> 0.5 & rnorm>fnorm) - warning("zgscal(tzero): GCG iteration failed; solving with pinv"); + if (xnorm> 0.5 && rnorm > fnorm) + warning ("zgscal(tzero): GCG iteration failed; solving with pinv"); ## perform brute force least squares; construct F - Am = eye(nmp); - for ii=1:nmp - Am(:,ii) = zgfmul(a,b,c,d,Am(:,ii)); + Am = eye (nmp); + for ii = 1:nmp + Am(:,ii) = zgfmul (a, b, c, d, Am(:,ii)); endfor ## now solve with qr factorization - x = pinv(Am)*z; + x = pinv (Am) * z; endif + endfunction