From octave-sources-request at bevo dot che dot wisc dot edu Tue Mar 9 13:01:44 1999 Subject: Various special matrices. From: peda at sectra dot se To: octave-sources at bevo dot che dot wisc dot edu Date: Tue, 9 Mar 1999 20:25:26 +0100 Hi! I took some time to implement the following functions found in Matlab: wilkinson, hadamard, pascal, rosser and compan As far as I can see, they do the same thing as their Matlab counterpart, except for hadamard. My version handles fewer cases than the Matlab version. /Peter -----------------8<-------------------------- ## Copyright (C) 1999 Peter Ekberg ## ## This file is not yet part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 2, or (at your option) ## any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA ## 02111-1307, USA. ## usage: wilkinson (n) ## ## Return the Wilkinson matrix of order n. ## ## See also: hankel, vander, sylvester_matrix, hilb, invhilb, toeplitz ## hadamard, rosser, compan, pascal ## Author: peda function retval = wilkinson (n) if (nargin != 1) usage ("wilkinson (n)"); endif nmax = length (n); if ~(nmax == 1) error ("wilkinson: expecting scalar argument, found something else"); endif side = ones(n-1,1); center = abs(-(n-1)/2:(n-1)/2); retval = diag(side, -1) + diag(center) + diag(side, 1); endfunction -----------------8<-------------------------- -----------------8<-------------------------- ## Copyright (C) 1999 Peter Ekberg ## ## This file is not yet part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 2, or (at your option) ## any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA ## 02111-1307, USA. ## usage: hadamard (n) ## ## Return the Hadamard matrix of order n. A matrix H is a Hadamard ## matrix if its elements are -1 or 1 and if H'*H = n*eye(n). ## This function only handles the cases where n is a power of 2. ## Hadamard matrices exist for n>2 where rem(n,4) = 0 and for n=1 ## and n=2. ## ## See also: hankel, vander, sylvester_matrix, hilb, invhilb, toeplitz ## wilkinson, rosser, compan, pascal ## Author: peda function retval = hadamard (n) if (nargin != 1) usage ("hadamard (n)"); endif nmax = length (n); if ~(nmax == 1) error ("hadamard: expecting scalar argument, found something else"); endif if ~(log2(n) == nextpow2(n)) error ("hadamard: expecting power of 2 argument, found something else"); endif if (n == 1) retval = 1; else part = hadamard(n/2); retval = [part part; part -part]; endif endfunction -----------------8<-------------------------- -----------------8<-------------------------- ## Copyright (C) 1999 Peter Ekberg ## ## This file is not yet part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 2, or (at your option) ## any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA ## 02111-1307, USA. ## usage: pascal (n, t) ## ## Return the Pascal matrix of order n if t=0. ## t defaults to 0. ## Return lower triangular Cholesky factor of the Pascal matrix if t=1. ## Return a transposed and permuted version of pascal(n,1) if t=2. ## ## pascal(n,1)^2 == eye(n) ## pascal(n,2)^3 == eye(n) ## ## See also: hankel, vander, sylvester_matrix, hilb, invhilb, toeplitz ## hadamard, wilkinson, rosser, compan ## Author: peda function retval = pascal (n, t) if (nargin > 2) || (nargin == 0) usage ("pascal (n, t)"); endif if (nargin == 1) t = 0; endif if !is_scalar (n) || !is_scalar (t) error ("compan: expecting scalar arguments, found something else"); endif retval = diag((-1).^[0:n-1]); retval(:,1) = ones(n, 1); for j=2:n-1 for i=j+1:n retval(i,j) = retval(i-1,j) - retval(i-1,j-1); endfor endfor if (t==0) retval = retval*retval'; elseif (t==2) retval = retval'; retval = retval(n:-1:1,:); retval(:,n) = -retval(:,n); retval(n,:) = -retval(n,:); if (rem(n,2) != 1) retval = -retval; endif endif endfunction -----------------8<-------------------------- -----------------8<-------------------------- ## Copyright (C) 1999 Peter Ekberg ## ## This file is not yet part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 2, or (at your option) ## any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA ## 02111-1307, USA. ## usage: rosser ## ## Return the Rosser matrix. ## ## See also: hankel, vander, sylvester_matrix, hilb, invhilb, toeplitz ## hadamard, wilkinson, compan, pascal ## Author: peda function retval = rosser if (nargin != 0) usage ("rosser"); endif retval = [ 611 196 -192 407 -8 -52 -49 29; 196 899 113 -192 -71 -43 -8 -44; -192 113 899 196 61 49 8 52; 407 -192 196 611 8 44 59 -23; -8 -71 61 8 411 -599 208 208; -52 -43 49 44 -599 411 208 208; -49 -8 8 59 208 208 99 -911; 29 -44 52 -23 208 208 -911 99; ]; endfunction -----------------8<-------------------------- -----------------8<-------------------------- ## Copyright (C) 1999 Peter Ekberg ## ## This file is not yet part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 2, or (at your option) ## any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA ## 02111-1307, USA. ## usage: compan (p) ## ## Return the companion matrix to the polynomial with coefficients p. ## ## See also: hankel, vander, sylvester_matrix, hilb, invhilb, toeplitz ## hadamard, wilkinson, rosser, pascal ## Author: peda function retval = compan (p) if (nargin != 1) usage ("compan (p)"); endif if !is_scalar (p) if !is_vector (p) error ("compan: expecting vector argument, found something else"); endif n = length (p); retval = diag(ones(n-2, 1), -1); retval(1,1:n-1) = -reshape(p(2:n), 1, n-1)/p(1); else retval = zeros(0); endif endfunction -----------------8<-------------------------- ________________________________ Peter Ekberg Sectra Communications AB - DBS Phone: +46 (0)35 - 17 84 94 GSM: +46 (0)70 - 37 886 73 Mail: peda at sectra dot se