toeblitz-1.0/0000755000442500044270000000000012206220510014242 5ustar cunninghamcunninghamtoeblitz-1.0/results/0000775000442500044270000000000012202602046015751 5ustar cunninghamcunninghamtoeblitz-1.0/startup.m0000775000442500044270000000340012202534277016142 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This is a startup file which initializes Matlab's settings to allow for better use of the toeplitz package. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% path(sprintf('%s/util', pwd), path) % Adds the subdirectories of the root folder to the path, allowing us to call functions from them. path(sprintf('%s/test', pwd), path) path(sprintf('%s/decomp', pwd), path) path(sprintf('%s/solve', pwd), path) path(sprintf('%s/trace', pwd), path) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The following code checks for the relevant MEX files (such as .mexa64 % or .mexglx, depending on the machine architecture), and it creates the % mex file if it can not find the right one. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Create the mex file if necessary. startup_checkMEX('decomp', 'decompToeplitzFastZohar') % Toeplitz inversion startup_checkMEX('decomp', 'logdetToeplitzFastZohar') % Toeplitz log determinant startup_checkMEX('trace', 'traceToeplitzFast') % Toeplitz trace %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %This simply clears all variables and closes all windows opened by Matlab. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all % Clears all variables. close all % Closes all windows opened by Matlab.toeblitz-1.0/trace/0000775000442500044270000000000012206220476015355 5ustar cunninghamcunninghamtoeblitz-1.0/trace/traceToeplitzFast.mexmaci640000775000442500044270000002040012202534277022545 0ustar cunninghamcunningham__TEXT__text__TEXT  __stubs__TEXT((__stub_helper__TEXT<$<__unwind_info__TEXT`P`__eh_frame__TEXTP__DATA__dyld__DATA__la_symbol_ptr__DATAH__LINKEDIT    h P  ce4X@d$  0@rpath/libmx.dylib 0@rpath/libmex.dylib 0@rpath/libmat.dylib 84/usr/lib/libstdc++.6.dylib 8/usr/lib/libgcc_s.1.dylib 8{/usr/lib/libSystem.B.dylib& ) ASLWAS%O%NUHAWAVAUATSHHIH;OIf=IH;>IH{2IH{&HI>AEGDAAw?fAщHcAXAuMcCYDXD)At fA~zAWA;ffAȉHcAXAuHcYXt0ff1AOf.HcAXAu EfAYEXH[A\A]A^A_]%%%LXLLL@ 44)4  XzRx 4h  __<HT    . ;QY dyld_stub_binding_helper__dyld_func_lookup_mexFunction_mxCreateDoubleScalar_mxGetN_mxGetPrtoeblitz-1.0/trace/traceToeplitzFast.mex0000775000442500044270000002103012202534277021541 0ustar cunninghamcunninghamELF44 (# HH$$QtdRtdttGNURXRt0ο 'ұSa     `  |CEqX4 + "gKRm Z a ?z l x __gmon_start___init_fini__cxa_finalize_Jv_RegisterClassesmexFunctionmxGetNmxCreateDoubleScalarmxGetPrliboctinterp.soliboctave.solibcruft.soliblapack.so.3gflibblas.so.3gflibfftw3.so.3libfftw3f.so.3libreadline.so.6libncurses.so.5libdl.so.2libhdf5.so.6libz.so.1libgfortran.so.3libstdc++.so.6libm.so.6libgcc_s.so.1libc.so.6_edata__bss_start_endGLIBC_2.1.3Psi r      US[|tX[ hhhhh UVS$u]t $()9s ((9rƃ$[^]US.ktt $Ѓ[]Ë$ÐUWVS\u!]} $$E]$EԋF$E̋F$Eċ$}EEEE}̃EЋEMЃU׋UEډU܋Uԍ‰E䍶U19|E܃E9uuًUԾ}U}؉UU19|m9uuًU1Mut&9uŰE \[^_]ËE}EِUVSmêt&Ћu[^]US[ppY[o)8BP l xo`l ~ (D$ oooo GCC: (Ubuntu 4.4.3-4ubuntu5.1) 4.4.3ztmexFunctionp!: int5upk:P(P,(\y AmuT Tr6\ muL muD 6[ amuH n( i( j( V = b h = 6% $ > $ > : ; $ > .? : ; ' @: ; I : ; I 4: ; I 4: ; I 4: ; I  I &I /home/cunni/toeplitz/wemmick/toeplitz/trace/usr/include/octave-3.2.3/octavetraceToeplitzFast.cmex.h7AU?עp@mmf>>,jy.y.mJb.>,kx.x.b-<>,kx | zCB Ilong long intunsigned charprhs/home/cunni/toeplitz/wemmick/toeplitz/trace/traceToeplitzFast.cTr_Middlelong long unsigned intmxArrayGNU C 4.4.3nlhsT_colT_rowplhsshort unsigned intlong doublenrhsshort intmexFunctiontemp_diag_sumttzuW2u 24W4gu gzW V2u24V4gugzVAdV2V4AVP2P4YP.symtab.strtab.shstrtab.note.gnu.build-id.gnu.hash.dynsym.dynstr.gnu.version.gnu.version_r.rel.dyn.rel.plt.init.text.fini.eh_frame.ctors.dtors.jcr.dynamic.got.got.plt.data.bss.comment.debug_aranges.debug_pubnames.debug_info.debug_abbrev.debug_line.debug_frame.debug_str.debug_loc$2H.o``<8 @ll~HoUo d $$ m DD( vll0q`|xxxH    0%= ]"t 4"0- 8"3 p `l$D l   x    (5 K Z h t@   x  ' , ; B V"rz ~l crtstuff.c__CTOR_LIST____DTOR_LIST____JCR_LIST____do_global_dtors_auxcompleted.7021dtor_idx.7023frame_dummy__CTOR_END____FRAME_END____JCR_END____do_global_ctors_auxtraceToeplitzFast.c__i686.get_pc_thunk.bx_GLOBAL_OFFSET_TABLE__DYNAMIC__DTOR_END____dso_handle_fini__bss_start_end__gmon_start___edata_Jv_RegisterClasses__cxa_finalize@@GLIBC_2.1.3mexFunctionmxGetPr_initmxGetNmxCreateDoubleScalartoeblitz-1.0/trace/traceToeplitzBackup.m0000775000442500044270000000315212202534277021521 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % % This generates the trace of the product of a Toeplitz matrix and any % arbitrary square matrix of compatible size. This is a backup function % in case the C/MEX implementation fails for any reason. This function % computes the trace from the contributions from different parts of the % original matrices. In this case both the first row and column vectors % have been specified. % % The calling syntax is % % traceToeplitzFast(A, T_row, T_col) % % where % A is the arbitrary matrix, % T_row is the first row of the Toeplitz matrix, and % T_col is the first column of the Toeplitz matrix, written as a row vector. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ Tr ] = traceToeplitzBackup(A, T_row, T_col) Tr = 0; % Initializes the trace at zero. for i = 2:size(T_col, 2), Tr = Tr + T_col(1, i)*sum(diag(A, i-1)); % Adds the contributions from each diagonal of the bottom half of the toeplitz matrix. Tr = Tr + T_row(1, i)*sum(diag(A, -i+1)); % Adds the contributions from each diagonal of the upper triangle of the toeplitz matrix. end Tr = Tr + T_row(1, 1)*sum(diag(A)); % Adds the contributions from the main diagonal of the toeplitz matrix. endtoeblitz-1.0/trace/traceToeplitzFast.c0000775000442500044270000000467512202534277021212 0ustar cunninghamcunningham/*================================================================= * Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices * Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) *=================================================================*/ /* William Zhang, 2013 * This simple C function simply computes the trace of the product of a Toeplitz matrix T and an arbitrary matrix A. * Mathematically, it is identical to traceToeplitzBackup.m. The calling syntax is * * traceToeplitzFast(A, T_row, T_col) * * where * A is the arbitrary matrix, * T_row is the first row of the Toeplitz matrix, and * T_col is the first column of the Toeplitz matrix, written as a row vector. */ #include #include /* Input Arguments */ #define A_IN prhs[0] #define T_row_IN prhs[1] #define T_col_IN prhs[2] /* Output Arguments */ #define Tr_OUT plhs[0] void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) { double *A, Tr, *T_row, *T_col, temp_diag_sum, *Tr_Middle; int n, i, j; /*Get the size of the matrices we're dealing with.*/ n = mxGetN(A_IN); /* Create a scalar for the return argument */ Tr_OUT = mxCreateDoubleScalar(0); /*Assign pointers to variables we will use to pass along inputs and outputs.*/ A = mxGetPr(A_IN); T_row = mxGetPr(T_row_IN); T_col = mxGetPr(T_col_IN); Tr_Middle = mxGetPr(Tr_OUT); /*Add the contribution from the top half of the Toeplitz matrix (and therefore the bottom half of the arbitrary matrix A), not counting the principal diagonal.*/ for (i = 1; i < n; i++) { temp_diag_sum = 0; for (j = 0; j < i; j++) { temp_diag_sum += A[(n-1)*n - (i-1)*n + j*(n+1)]; } Tr += temp_diag_sum*T_row[n-i]; } /*Add the contribution from the bottom half of the Toeplitz matrix (and therefore the top half of the arbitrary matrix A), not counting the principal diagonal.*/ for (i = 1; i < n; i++) { temp_diag_sum = 0; for (j = 0; j < i; j++) { temp_diag_sum += A[n-1 - (i-1) + j*(n+1)]; } Tr += temp_diag_sum*T_col[n-i]; } /*Add the contribution from the principal diagonal of the Toeplitz matrix defined by T_row and coL_vec, and the principal diagonal of arbitrary matrix A.*/ for (i = 0; i < 1; i++) { temp_diag_sum = 0; for (j = 0; j < n; j++) { temp_diag_sum += A[j*(n+1)]; } Tr += T_row[0]*temp_diag_sum; *Tr_Middle = Tr; } return; } toeblitz-1.0/trace/traceToeplitzFast.mexa640000775000442500044270000001507512202534277022070 0ustar cunninghamcunninghamELF>@@@ @8@ll pp p    $$PtdQtdGNUh.%ca:UCOo@  5 [?  F"3 __gmon_start____cxa_finalize_Jv_RegisterClassesmexFunctionmxGetNmxCreateDoubleScalarmxGetPrlibmx.solibmex.solibmat.solibm.so.6libstdc++.so.6libpthread.so.0libc.so.6traceToeplitzFast.mexa64MEXGLIBC_2.2.5DB4Qui         H_H5 % @% h% h% h% hHHU HtHÐU= HATSubH=8 t H= H L% H] L)HHH9s DHH= AH2 H9r [A\fH= UHtH HtH= @ÐUHAWAVAUATSH8HuIH9IAfWHUHI?HIIIIHEH8HEAEAAT$UAEžA|$AA3HcX˃9uD)HAYXEEȃEA9~LME~DDfWDEA|$AAM0HcX˃9uIcAYXEEȃAA9~LME~DDLUEfWE~%At$fWHcX˃A9AYXEHEH8[A\A]A^A_ÐUHSHH HtH HHHuH[ÐH_H;0zRx ,AC T  clv  o  `x` oooooo &6GCC: (GNU) 4.4.6 20120305 (Red Hat 4.4.6-4)GCC: (GNU) 4.4.7 20120313 (Red Hat 4.4.7-3).symtab.strtab.shstrtab.note.gnu.build-id.gnu.hash.dynsym.dynstr.gnu.version.gnu.version_d.gnu.version_r.rela.dyn.rela.plt.init.text.fini.eh_frame_hdr.eh_frame.ctors.dtors.jcr.data.rel.ro.dynamic.got.got.plt.bss.comment$.o(8 @HoUo8do s`}xx` P@@  Lp p      8 0 XH @X1  x   @  p         @p * 8 E `[ j x x h         # , 2:A P dh}" call_gmon_startcrtstuff.c__CTOR_LIST____DTOR_LIST____JCR_LIST____do_global_dtors_auxcompleted.6349dtor_idx.6351frame_dummy__CTOR_END____FRAME_END____JCR_END____do_global_ctors_auxtraceToeplitzFast.c_fini_GLOBAL_OFFSET_TABLE___dso_handle__DTOR_END____bss_start_end_edata_DYNAMIC_initmxGetPrmxGetN__gmon_start___Jv_RegisterClassesMEXmxCreateDoubleScalar__cxa_finalize@@GLIBC_2.2.5mexFunctiontoeblitz-1.0/trace/traceToeplitzFast.mexw640000775000442500044270000001700012202534277022104 0ustar cunninghamcunninghamMZ@ !L!This program cannot be run in DOS mode. $3Gw)Lw)Lw)LLu)LLu)L~Lr)Lw(Ln)LL})LLv)LLv)LRichw)LPEd:| R"  xpq@0%W$"PP@` .textb   `.rdata @@.data`0@.pdata@@@.rsrcP@@.reloc `@BH\$UVWATAUAVAWHPI )t$@IHQHfWf(;HH*HOLH$HOLHD$8HLHD$($AHcE3LHD$0EL;NH$LLIBIIH+LHjHHILHHILOd"H$HHL$ fDf(I|vIHZJL)I+HMHxHpIAHHHHLfffffAXLXDXXDHHuH$H$M;}%HBMBIHH II+XIHuAY$HL$ IIMHIAXL;(L$Lt$(L|$0ESHL;H<MJt7M+f(II|UHRLJLI+IL+ILIBM4HHHAXMBXXXD9IHuE3I;})HBLGHI+HILIH+XIHuYIHXL;XH|THBHJLJHLHIII+HMLAX0MX4(X0BXtIHuL;}$HBLII+ILX1IHuHD$8H$Y0XA7(t$@HPA_A^A]A\_^]%^ %P %B @SH  HH HHHuC#H#H H 3H [HHXHhHxL` AUAVAWH 3ML8 #DeH%0HXH;t 3HuAt_H  LHH  MLHHI;rZH9}t H9EtHMh H HEH rL H ]H< L;uL;tLLI s H,H-=EH=3eH%0HXH;t 3Hut 7>H` H I uH' H  au HH[H9=lt!H ctMĺII H\$@Hl$HH|$PLd$XH A_A^A]HHXHpHxATH0ILXu9u 3ۉXtu7Hz HtЋ؉D$ LƋI0؉D$ LƋI؉D$ u5u1L3IL3IL Mt L3IAӅtu7LƋI#ˋىL$ tH HtLƋIЋ؉D$ 3ۉ\$ H\$@Ht$HH|$PH0A\H\$Ht$WH IHukLNjHH\$0Ht$8H _H i@SH HH |VHD$8Hu H~H N(HD$8H 4HD$@HHLD$@HT$8HHL$8HHL$@HSHH [H(GHH(H\$WH H{H=tHHtHH;rH\$0H _H\$WH HSH=LHHtHH;rH\$0H _HMZf9t3HcH PAPADDINGXXPADDINGPADDINGXXPADDINGPADDINGXXPADDINGPADDINGXXPADDINGPADDINGXXPAD toeblitz-1.0/logdetToeplitz.m0000775000442500044270000000265712202534277017466 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This is simply a wrapper function that calls decompToeplitz in a way that computes only the log determinant, and does not waste computations and memory constructing and storing the inverse. We can specify a runMode here to run decompToeplitz in the desired mode. Here are the accepted values of runMode: % runMode = 0 % This calls the C/MEX implementation of Zohar's modification of the Trench algorithm. % runMode = 1 % This calls the backup Matlab implementation of Zohar's modification of the Trench algorithm. % runMode = -1 % This calls the generic MATLAB function logdet. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[ld] = logdetToeplitz(T_row, runMode) if nargin < 2 % Handles the case when runMode is not specified. ld = decompToeplitz(T_row, 1); % This allows decompToeplitz to automatically choose the best method. else % Handles the case when runMode is specified. ld = decompToeplitz(T_row, 1, runMode); % This forces decompToeplitz to use the specified method. end endtoeblitz-1.0/solveToeplitz.m0000775000442500044270000000653412202534277017336 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This solves a system Tx = b for x, where x and b are column vectors and T is a toeplitz matrix. It uses the method specified by runMode, or chooses a method based on the size of the matrix and vector. % Accepted values of runMode are: % % runMode = 'LD' This directs the code to execute the Levinson-Durbin algorithm, which works in % runMode = 'PCG' This directs the code to execute the preconditioned conjugate gradient method with automatic choice of preconditioner. The code will always choose Chan. % runMode = 'Mldivide' This directs the code to use the built-in MATLAB mldivide method. % runMode = 'PCG_Chan' This explicitly directs the code to use the preconditioned conjugate gradient method with the Chan preconditioner. % runMode = 'PCG_Circulant_Embedding' This explicitly directs the code to use the preconditioned conjugate gradient method with a preconditioner produced by embedding our Toeplitz matrix in a circulant matrix and adding it to the other block. This is described more fully in the "Conjugate gradient methods for Toeplitz systems" 1996 article we reference in the paper accompanying this software package. % runMode = 'PCG_None' This explicitly directs the code to use the unconditioned conjugate gradient method. % % The calling syntax for solveToeplitz is % % solveToeplitz(T_row, b [, runMode]) % % where % % T_row is the first row of the Toeplitz matrix T % b is a the column vector % runMode (optional) is the method used to solve for x. % % The output of solveToeplitz is a single variable % % x, which is the solution to the matrix equation Tx = b. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x] = solveToeplitz(T_row, b, runMode) size_T_row = size(T_row, 2); % Finds the size of the Toeplitz matrix. This is used for determining which method to use (different methods are fastest at different matrix sizes). if nargin < 3 if size_T_row < 500 runMode = 'mldivide'; % Uses Mldivide for smaller matrices. elseif size_T_row >= 500 runMode = 'pcg'; % Uses the Chan-preconditioned conjugate gradient method for larger matrices. end end x = solveToeplitzMode(T_row, b, runMode); end function [x] = solveToeplitzMode(T_row, b, runMode) %This function links the runMode handles to the actual function which implement the different methods. if strcmp(lower(runMode), 'ld') x = solveToeplitzLD(T_row, b); elseif strcmp(lower(runMode), 'mldivide') x = solveToeplitzMldivide(T_row, b); elseif strcmp(lower(runMode(1, 1:3)), 'pcg') if length(runMode) == 3 %This case allows our code to choose the preconditioner automatically. x = solveToeplitzPCG(T_row, b, 'auto'); elseif runMode(1, 4) == '_' %This case forces solveToeplitzPCG to run with a specified preconditioner. x = solveToeplitzPCG(T_row, b, runMode(5:end)); else error('Invalid runMode value given.') end else error('Invalid runMode value given.') end endtoeblitz-1.0/README0000775000442500044270000003002112202603643015133 0ustar cunninghamcunningham==================================================================== William Zhang, Washington University in St. Louis John P. Cunningham, Columbia University Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices Copyright (C) 2013 William B. Zhang and John P. Cunningham This program 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 3 of the License, or (at your option) any later version. This program 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 this program. If not, see . ==================================================================== ==================================================================== Basic Usage Example ==================================================================== Here, we present a basic example of how to use our Toeblitz toolkit, with some random matrices and variables: T_row = genPDToeplitz(50); % This creates the variable T_row, the first row of a random 50 x 50 symmetric positive definite Toeplitz matrix T, which we will study. A = randn(50, 50); % This creates a random 50 x 50 matrix to use in computing the trace of a product of a Toeplitz matrix with an arbitrary matrix. b = randn(50, 1); % This creates a random 50 x 1 vector to use in the matrix Tx = b for testing solveToeplitz. tr = traceToeplitz(A, T_row); % This computes the trace of the product of A and T in a fast way. [invT, ldT] = invToeplitz(T_row); % This computes the inverse of T and the log determinant of T in a fast way. x = solveToeplitz(T_row, b); % This solves the system Tx = b for x. ==================================================================== ==================================================================== C/MEX: MEX functions are programs written in C, C++, or FORTRAN that are callable from Matlab. This allows bottlenecks in native Matlab code to be optimized through the use of a compiled language. In this package, we have pre-compiled Matlab MEX functions from C code for several popular operating systems. If we do not have a precompiled version for your operating system, the startup.m file will automatically attempt create one for your operating system. This requires a C compiler. Supported compilers can be found at . More information on C/MEX can be found at . Note: We have also pre-compiled MEX files for Octave, but have only done so (and tested our C functions) on an Ubuntu distribution of Linux. ==================================================================== ==================================================================== Startup.m This script will run by default if Matlab starts in the /toeplitz directory. If Matlab does not start in the /toeplitz directory, you should switch to it using "cd .../toeplitz", where .../toeplitz is the full directory path of the directory on your system. Then, you should run startup.m manually by simply typing "startup" into Matlab. It adds the locations of all the functions contained in our script to Matlab's search path, and attempts to create the relevant MEX files for your operating system, if it does not already exist. ==================================================================== ==================================================================== solveToeplitz.m solveToeplitz allows us to solve a Toeplitz system Tx = b for the unknown vector x. We call solveToeplitz with the following syntax: x = solveToeplitz(T_row, b [, runMode]) where b is the column vector b from our equation, T_row is the first row of the Toeplitz matrix T, x is the desired vector from our equation, and runMode is an optional argument forcing the code to use a particular method to obtain the solution x. Valid values of runMode are presented in the documentation in the main solveToeplitz function. ==================================================================== ==================================================================== logdetToeplitz.m logdetToeplitz allows us to compute the log determinant of a Toeplitz matrix T. logdetToeplitz is faster than obtaining the log determinant as a secondary output from invToeplitz, since it wastes no memory constructing the full inverse. It is linear in memory and quadratic in runtime. The calling syntax for logdetToeplitz is ldT = logdetToeplitz(T_row [, runMode]) where ldT is the log determinant of T, T_row is the first row of the Toeplitz matrix T, and runMode is an optional argument forcing the code to use a particular method to obtain the log determinant ldT. Valid values of runMode are presented in the documentation in the main logdetToeplitz function. ==================================================================== ==================================================================== traceToeplitz.m This function takes the trace of the product of a Toeplitz matrix T and an arbitrary matrix A in a smart way. The calling syntax is tr = traceToeplitz(A, T_row [, T_col, runMode]) where A is the arbitrary matrix, T_row is the first row of T, tr is the trace of the product A*T, T_col is an optional argument giving the first column of T (when T is asymmetric), and runMode is an optional argument forcing the code to use a particular method to obtain the trace. Valid values of runMode are presented in the documentation in the main traceToeplitz function. Note that in our tests, trace(A*T) is used as the MATLAB built-in function. sum(sum(transpose(A).*T)) is actually equivalent, and gives the answer faster than trace(A*T). However, our implementation exploits Toeplitz structure to perform the computation in the minimum (n^2 + n - 1) floating point operations, whereas sum(sum(transpose(A).*T)) requires (2n^2 - 1) operations. In our analysis we choose trace(A*T) as the comparison in order to illustrate the contrast between its O(n^3) complexity and the O(n^2) complexity of our traceToeplitz method. ==================================================================== ==================================================================== invToeplitz.m invToeplitz computes the inverse of a Toeplitz matrix T, and gives the log determinant of T for free in an optional second output argument. The calling syntax for invToeplitz is [invT [, ldT]] = invToeplitz(T_row [, runMode]) where invT is the inverse of T, ldT is an optional output argument giving the log determinant of T, T_row is the first row of the Toeplitz matrix T, and runMode is an optional argument forcing the code to use a particular method to obtain the inverse invT. Valid values of runMode are presented in the documentation in the main invToeplitz function. ==================================================================== ==================================================================== testToeplitz.m You can use this script to test the accuracy and speed of the supported functions in this package. The calling syntax is testToeplitz(method, sizes_vec, runs, save_flag)". where method is a string giving the name of the method to be tested (valid values of method are described in the documentation of the main function, testToeplitz), sizes_vec is a row vector of matrix sizes for which to test that method, runs is the number of iterations to average over for each data point, and save_flag is a Boolean indicating whether the results of the test should be saved, or simply displayed. To display the results without saving, save_flag can be set to 0 (false), or simply omitted. To save your results, save_flag must be set to 1. In this case, the figures will not be displayed, but will be saved as both .pdf files and .fig files to /toeplitz/results/XXX/, where XXX is a number assigned to that set of data, ranging from 001 to 999. ==================================================================== ==================================================================== recreateToeblitzPlots.m This script, when run with no input parameters, will generate and save to disk the three plots presented in the Toeblitz paper. It should be able to complete running overnight on a modern laptop. For a quicker verification, the -1 argument can be used, as follows: recreateToeblitzPlots(-1) Other options exist to only produce one of the three plots in the Toeblitz paper, and these are documented in the recreateToeblitzPlots.m file itself. ==================================================================== ==================================================================== convertToeplitz.m convertToeplitz is provided in the case that the user wishes to convert between our abbreviated representation of a Toeplitz matrix as its first row (and sometimes also its first column, in the asymmetric case) and the full explicit representation. Its calling syntax is T = convertToeplitz(T_row [, T_col]) [T_row [, T_col]] = convertToeplitz(T) where T is the explicit full representation of the Toeplitz matrix, and T_row and T_col are the first row and column, respectively. We see that T_col is an optional input used to specify asymmetric Toeplitz matrices, and that it is only given as an output when an asymmetric Toeplitz matrix is given as the input argument. ==================================================================== ==================================================================== The BibTeX citations for the primary papers used in this project are below: @book{golub1996matrix, title={Matrix computations}, author={Golub, Gene H and Van Loan, Charles F}, volume={3}, year={2012}, publisher={The Johns Hopkins University Press} } @book{rasmussen2006gaussian, title={Gaussian processes for machine learning}, author={Rasmussen, Carl E and Williams, Christopher KI}, year={2006}, publisher={MIT Press} } @article{levinson1947wiener, title={The {W}iener {R}{M}{S} (root mean square) error criterion in filter design and prediction}, author={Levinson, Norman}, journal = {Journal of Mathematics and Physics}, volume={25}, number={4}, pages={261–-278}, year={1947} publisher={MIT Press} } @article{durbin1960fitting, title={The fitting of time-series models}, author={Durbin, James}, journal={Revue de l'Institut International de Statistique}, pages={233--244}, year={1960}, publisher={JSTOR} } @article{holmes1989random, title={On random correlation matrices {I}{I}: the {T}oeplitz case}, author={Holmes, Robert B}, journal={Communications in Statistics-Simulation and Computation}, volume={18}, number={4}, pages={1511--1537}, year={1989}, publisher={Taylor \& Francis} } @article{zohar1969toeplitz, title={Toeplitz matrix inversion: the algorithm of {W}{F} {T}rench}, author={Zohar, Shalhav}, journal={Journal of the Association for Computing Machinery (JACM)}, volume={16}, number={4}, pages={592--601}, year={1969}, publisher={ACM} } @article{chan1996conjugate, title={Conjugate gradient methods for {T}oeplitz systems}, author={Chan, Raymond H and Ng, Michael K}, journal={SIAM review}, volume={38}, number={3}, pages={427--482}, year={1996}, publisher={SIAM} } @article{shewchuk1994introduction, title={An introduction to the conjugate gradient method without the agonizing pain}, author={Shewchuk, Jonathan R}, year={1994}, journal={CMU Technical Report} publisher={Carnegie Mellon University, Pittsburgh, PA} } @article{minka2003lightspeed, title={The lightspeed {MATLAB} toolbox}, author={Minka, Tom}, journal={Efficient operations for Matlab programming, Version}, volume={2}, year={2003} } ====================================================================toeblitz-1.0/util/0000775000442500044270000000000012202534277015237 5ustar cunninghamcunninghamtoeblitz-1.0/util/multiplyToeplitz.m0000775000442500044270000000307312202534277021035 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This uses a fast Fourier transform to compute the product y = Tx, where T is a toeplitz matrix and x is a column vector of appropriate size. T is embedded in a circulant matrix and manipulated according to page 202 of "Matrix Computations", 3rd edition by Gene H. Golub. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [y] = multiplyToeplitz(T_row, x) n = size(T_row, 2); % Gets the dimension of the Toeplitz matrix. embed_circulant = T_row(1, end:-1:2); % Extracts information from first row of Toeplitz matrix. embed_circulant = vertcat(transpose(T_row(1, 1:end)), embed_circulant(:)); % Combines information from first column of Toeplitz matrix with the information from the first row extracted in the previous step. Note that these two steps can be condensed into a single faster step if the toeplitz matrix is known to be symmetric. y = ifft(fft(embed_circulant).*fft(vertcat(x, zeros(n-1, 1)))); % Computes the result of the circulant matrix multiplication in a fast way. y = y(1:n, 1); % Discards the extra trailing elements in the circulant matrix result, giving us the desired product of a Toeplitz matrix and a vector. endtoeblitz-1.0/util/genPDToeplitz.m0000775000442500044270000000520212202534277020147 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This function generates the first row of a random symmetric positive definite Toeplitz matrix of size n. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This is the main function. It generates a random positive definite Toeplitz matrix of size n from the algorithm given on page 4 of "On random correlation matrices: II. The Toeplitz case" (1989) by R. B. Holmes at the Lincoln Laboratories at the Massachusetts Institute of Technology. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [T_row] = genPDToeplitz(n) % Generate a random vector from the uniform distribution on the unit n-sphere. rand_vec = randn(n, 1); % Generate a random n-dimensional vector. This essentially generates the direction of our vector. rand_vec = rand_vec/norm(rand_vec, 2); % Normalize our direction to the desired magnitude, 1. % Now generate a positive definite Toeplitz matrix as a deterministic function of our psuedorandom vector. for j = 1:n % Loop over the entries in the first row of the Toeplitz matrix. T_row(1, j) = phi_k(abs(1-j), n, rand_vec); % Produce the first row of the Toeplitz matrix. end rand_mat_mag = abs(randn(1, 1)); % Produce a random positive number by which to scale this matrix (just for some additional variation). T_row = rand_mat_mag*T_row; % Scale the first row by this random positive number. end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This is a subfunction which generates an element of a positive definite Toeplitz matrix, given an integer k = abs(i-j), the absolute value of the difference between the indices of the entry in the matrix, an integer n, the desired size of the matrix, and rand_vec, a n x 1 random vector in the unit n-sphere. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[new_x] = phi_k(k, n, rand_vec) rand_vec(n+1:n+k+1, 1) = 0; % Extend our random vector with k more zeros. new_x = transpose(rand_vec(1:n, 1))*rand_vec(1+k:n+k, 1); % Compute the required sum in a vectorized way (from page 4 of "On random correlation matrices: II. The Toeplitz case" (1989) by R. B. Holmes). endtoeblitz-1.0/util/multiplyinvCirculant.m0000775000442500044270000000203012202534277021654 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This function solves the matrix equation C*x = b for x, where C is a circulant matrix. This is equivalent to left-multiplying b by the inverse of C, C^(-1). It takes c, the first column of c, as input. C is never explicitly required or stored. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x] = multiplyinvCirculant(c, b) x = ifft(fft(b)./fft(c)); % Performs the computation. The theoretical justification of this is a consequence of the fast Toeplitz multiplication method on page 202 of "Matrix Computations", 3rd edition by Gene H. Golub, and the properties of the Fourier transform. endtoeblitz-1.0/util/recreateToeblitzPlots.m0000775000442500044270000006742612202546227021766 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This function makes it convenient to reproduce all the plots included in the accompanying journal article at once. It's essentially a modified version of testToeplitz with no input arguments. You can just call it with no input variables and let it run. runMode is an optional variable, which if set to -1, will run the code for smaller matrices, making it easy to verify the code in a couple minutes, and allows one to run the code many times to find errors. If runMode is set to 0, it doesn't do anything, and the resulting behavior is equivalent to omitting it to begin with. Setting runMode to 1, 2, or 3 will run only a particular portion of the code needed to create one particular graph. A second optional input argument, fig_save, specifies whether we save the produced plots as MATLAB .fig files in addition to .pdf files. In order to specify the second input argument, runMode must be explicitly specified. fig_save defaults to 0 (only saves .pdf files) if not specified. % % In total, this function creates four plots when called in runMode 0 (default) or -1: one accuracy plot on a log-linear axis; one plot of iterations on a linear axis, and two loglog plots for complexity, one on a true loglog axis, and the other a linear plot of the logs of the two variables, used for regressions. When called in runMode 1, only the one accuracy plot is produced; when called in runMode 2, only the one iterations plot is produced; and when called in runMode 3, only the two complexity plots are produced. % % Calling syntax: % % recreateToeblitzPlots This simply runs the full recreateToeblitzPlots procedure. It should complete running overnight on a modern laptop computer. % recreateToeblitzPlots() This is equivalent to simply calling recreateToeblitzPlots. % recreateToeblitzPlots(0) This is equivalent to simply calling recreateToeblitzPlots. % recreateToeblitzPlots(-1) This runs a shorter version of the code to produce graphs for quick verification. % recreateToeblitzPlots(1) This produces only the accuracy check (full version) plot. % recreateToeblitzPlots(2) This produces only the iteration check (full version) for the various preconditioners. % recreateToeblitzPlots(3) This produces only the complexity testing (full version) plot for the various methods. % recreateToeblitzPlots(0, 0) This is equivalent to calling recreateToeblitzPlots(0) or simply recreateToeblitzPlots. % recreateToeblitzPlots(0, 1) This saves the produced figures as MATLAB .fig files in addition to .pdf files for editing and presentation. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [] = recreateToeblitzPlots(runMode, fig_save) % Set runMode to 0 by default. if nargin < 1 runMode = 0; end % Set fig_save to 0 by default. if nargin < 2 fig_save = 0; end % The following code sets default values for runs and save_flag if they are not specified. save_flag = 1; % Set save_flag to true. runs = 10; % Set the number of runs to 10. method = 'Placeholder'; % Set a placeholder for method at first. % For testing only: if runMode is -1 (quick verification), then only do a single run for each case, and do not save the output; instead, display it immediately. if runMode == -1 runs = 3; save_flag = 0; end % The following code finds a folder to save the result plots in, if save_flag is set to 1. folder = '0'; if save_flag for i = 1:999 % Searches from 001 to 999 for an unused folder name of the form 'test_XXX', where XXX is the three-digit number. if folder == '0' % Runs the following code if we have not yet found a suitable folder. test_folder = sprintf('./results/test_%03d/', i); % Finds the subdirectory in which the files will be placed. if not(exist(test_folder)) % Checks to see if there is already a folder with that name in the ./results/ directory. If not, it creates the folder and proceeds. folder = test_folder; % Changes the value of folder, stopping future checks for folder names. mkdir(folder); % Actually creates the directory. fprintf(['Printing to directory ' folder '.\n']) % Informs the user where the output will be saved. end end end end if or(or(runMode == 0, runMode == -1), runMode == 1) method = 'accuracyToeplitz'; % Sets the method to traceToeplitz. sizes_vec = [2:100:4502]; % Sets up a sizes_vec for small matrices. This will let us do a quick accuracy check. end % For testing only: if runMode == -1 sizes_vec = [100:20:200]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The following code tests the accuracy of our methods. Routines for testing other functions work similarly. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if strcmp(method, 'accuracyToeplitz') % This checks if 'accuracyToeplitz' is the method specified for testing. accuracy_strs = cell(4, 1); % This creates a cell array of the names of our different methods for use in labeling plots and legends. accuracy_strs(1, 1) = cellstr('solveToeplitz'); accuracy_strs(2, 1) = cellstr('logdetToeplitz'); accuracy_strs(3, 1) = cellstr('invToeplitz'); accuracy_strs(4, 1) = cellstr('traceToeplitz'); vals_times = zeros(4, length(sizes_vec)); % This initializes the matrix in which all our test results will be (temporarily) stored. for i = 1:length(sizes_vec) % This loops over each the sizes of matrices we want to test. current_size_str = num2str(sizes_vec(i)); % Finds the current matrix size being worked on. new_vals_times = zeros(4, 1); % This initializes new_vals_times (and resets it for each size of matrix tested). for j = 1:runs % For each size matrix, this loops 'runs' times to obtain an average result. n = sizes_vec(i); run_str = num2str(j); % Finds the current run being worked on. fprintf(['Now on run number ' run_str ' of matrix size ' current_size_str ' for method ' method '.\n']) % Informs the user where we are in the program. T_row = genPDToeplitz(n); % This produces the first row of a (more or less) random positive definite Toeplitz matrix. b = transpose([1:n]); A = randn(n, n); new_vals_times = new_vals_times + test_accuracyToeplitz(T_row, A, b); % This calculates the trace of the product with an arbitrary matrix in various ways, and adds the timings and accuracy metrics to a cumulative total in new_vals_times (see test_traceToeplitz documentation for more detail). end vals_times(1:4, i) = new_vals_times(1:4, 1)/runs; % This saves the results in new_vals_times to the appropriate column in vals_times, dividing by the number of runs to obtain an average. end % Show or save results: These are results for accuracy. test_plot('semilogy', sizes_vec, vals_times(1:1:4, 1:end), accuracy_strs, 'Norm difference vs. built-in MATLAB methods', save_flag, folder, 'test_accuracy', fig_save); fprintf(['Printing to directory ' folder '.\n']) % Informs the user where the output will be saved. end if or(or(runMode == 0, runMode == -1), runMode == 2) method = 'iterToeplitz'; sizes_vec = [500:100:4500]; end % For testing only: if runMode == -1 sizes_vec = [50:50:250]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The following code tests the number of iterations required for various preconditioned conjugate gradient methods. Routines for testing other functions work similarly. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if strcmp(method, 'iterToeplitz') % This checks if 'iterToeplitz' is the method specified for testing. iter_strs = cell(3, 1); % This creates a cell array of the names of our different methods for use in labeling plots and legends. iter_strs(1, 1) = cellstr('Unconditioned'); iter_strs(2, 1) = cellstr('Chan'); iter_strs(3, 1) = cellstr('Circulant Embedding'); vals_times = zeros(3, length(sizes_vec)); % This initializes the matrix in which all our test results will be (temporarily) stored. for i = 1:length(sizes_vec) % This loops over each the sizes of matrices we want to test. current_size_str = num2str(sizes_vec(i)); % Finds the current matrix size being worked on. new_vals_times = zeros(3, 1); % This initializes new_vals_times (and resets it for each size of matrix tested). for j = 1:runs % For each size matrix, this loops 'runs' times to obtain an average result. n = sizes_vec(i); run_str = num2str(j); % Finds the current run being worked on. fprintf(['Now on run number ' run_str ' of matrix size ' current_size_str ' for method ' method '.\n']) % Informs the user where we are in the program. T_row = genPDToeplitz(n); % This produces the first row of a (more or less) random positive definite Toeplitz matrix. b = transpose([1:n]); A = randn(n, n); new_vals_times = new_vals_times + test_iterToeplitz(T_row, b); % This calculates the trace of the product with an arbitrary matrix in various ways, and adds the timings and accuracy metrics to a cumulative total in new_vals_times (see test_traceToeplitz documentation for more detail). end vals_times(1:3, i) = new_vals_times(1:3, 1)/runs; % This saves the results in new_vals_times to the appropriate column in vals_times, dividing by the number of runs to obtain an average. end % Show or save results: These are results for accuracy. test_plot('linear', sizes_vec, vals_times(1:1:3, 1:end), iter_strs, 'Number of PCG Iterations', save_flag, folder, 'test_iter', fig_save); fprintf(['Printing to directory ' folder '.\n']) % Informs the user where the output will be saved. end if or(or(runMode == 0, runMode == -1), runMode == 3) method = 'complexityToeplitz'; % Sets the method to speedToeplitz. sizes_vec_small = floor(10.^[2.7:0.1:3.9]); % Sets up a sizes_vec_small for matrices large enough so that the computational time is measureable, but not so large that they break MATLAB. This will let us do a quick accuracy check. sizes_vec_large = floor(10.^[4:0.1:5]); % Sets up a sizes_vec_large for matrices large enough so that the computational time is measureable, but not so large that they break MATLAB. end % For testing only: if runMode == -1 sizes_vec_small = floor(10.^[2.7:0.2:3.5]); sizes_vec_large = floor(10.^[3.9:0.2:4.3]); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The following code tests the speed of our methods. Routines for testing other functions work similarly. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if strcmp(method, 'complexityToeplitz') % This checks if 'complexityToeplitz' is the method specified for testing. complexity_strs = cell(8, 1); % This creates a cell array of the names of our different methods for use in labeling plots and legends. complexity_strs(1, 1) = cellstr('solveToeplitz'); complexity_strs(2, 1) = cellstr('logdetToeplitz'); complexity_strs(3, 1) = cellstr('invToeplitz'); complexity_strs(4, 1) = cellstr('traceToeplitz'); complexity_strs(5, 1) = cellstr('MATLAB mldivide'); complexity_strs(6, 1) = cellstr('MATLAB logdet'); complexity_strs(7, 1) = cellstr('MATLAB inv'); complexity_strs(8, 1) = cellstr('MATLAB trace'); vals_times_small = zeros(8, length(sizes_vec_small)); % This initializes the matrix in which our test results for smaller matrices will be (temporarily) stored. vals_times_large = zeros(2, length(sizes_vec_large)); % This initializes the matrix in which our test results for larger matrices will be (temporarily) stored. for i = 1:length(sizes_vec_small) % This loops over each the sizes of matrices we want to test. n = sizes_vec_small(i); current_size_str = num2str(n); % Finds the current matrix size being worked on. new_vals_times = zeros(8, 1); % This initializes new_vals_times (and resets it for each size of matrix tested). for j = 1:runs % For each size matrix, this loops 'runs' times to obtain an average result. run_str = num2str(j); % Finds the current run being worked on. fprintf(['Now on run number ' run_str ' of matrix size ' current_size_str ' for method ' method '.\n']) % Informs the user where we are in the program. T_row = genPDToeplitz(n); % This produces the first row of a (more or less) random positive definite Toeplitz matrix. b = transpose([1:n]); A = randn(n, n); new_vals_times = new_vals_times + test_complexityToeplitz_small(T_row, A, b); % This calculates the trace of the product with an arbitrary matrix in various ways, and adds the timings and accuracy metrics to a cumulative total in new_vals_times (see test_traceToeplitz documentation for more detail). end vals_times_small(1:8, i) = new_vals_times(1:8, 1)/runs; % This saves the results in new_vals_times to the appropriate column in vals_times, dividing by the number of runs to obtain an average. end for i = 1:length(sizes_vec_large) % This loops over each the sizes of matrices we want to test. n = sizes_vec_large(i); current_size_str = num2str(n); % Finds the current matrix size being worked on. new_vals_times = zeros(2, 1); % This initializes new_vals_times (and resets it for each size of matrix tested). for j = 1:runs % For each size matrix, this loops 'runs' times to obtain an average result. run_str = num2str(j); % Finds the current run being worked on. fprintf(['Now on run number ' run_str ' of matrix size ' current_size_str ' for method ' method '.\n']) % Informs the user where we are in the program. T_row = genPDToeplitz(n); % This produces the first row of a (more or less) random positive definite Toeplitz matrix. b = transpose([1:n]); new_vals_times = new_vals_times + test_complexityToeplitz_large(T_row, b); % This calculates the trace of the product with an arbitrary matrix in various ways, and adds the timings and accuracy metrics to a cumulative total in new_vals_times (see test_traceToeplitz documentation for more detail). end vals_times_large(1:2, i) = new_vals_times(1:2, 1)/runs; % This saves the results in new_vals_times to the appropriate column in vals_times, dividing by the number of runs to obtain an average. end % Redistribute the data so that each actual data series is together in one row. vals_times_large = horzcat(vals_times_small(1:2, 1:end), vals_times_large); vals_times_small = vals_times_small(3:8, 1:end); % Show or save results: These are results for computational complexity. grand_loglogplot(sizes_vec_small, sizes_vec_large, vals_times_small, vals_times_large, complexity_strs, 'Computational Time', save_flag, folder, 'test_complexity', runMode, fig_save); grand_regressionplot(sizes_vec_small, sizes_vec_large, vals_times_small, vals_times_large, complexity_strs, 'Computational Time', save_flag, folder, 'test_complexityregression', runMode, fig_save); fprintf(['Printing to directory ' folder '.\n']) % Informs the user where the output will be saved. end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % test_accuracyToeplitz tests the accuracy of our Toeblitz method against built-in MATLAB functions and outputs the results in a column vector, new_vals_times. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[new_vals_times] = test_accuracyToeplitz(T_row, A, b) new_vals_times = zeros(4, 1); % Intialize new_vals_times as a vector of zeros. % Check the accuracy of our preconditioned conjugate gradient method. x = solveToeplitzMldivide(T_row, b); % Get a standard against which to measure the accuracy of our Chan preconditioned conjugate gradient method. xPCG = solveToeplitzPCG(T_row, b, 'Chan'); % Get the result using our method. new_vals_times(1, 1) = norm(x(:) - xPCG(:), 2); % Compares accuracy to built-in Matlab algorithm. % Check the accuracy of the inverse in our decompToeplitz method. [TiInvT, ldT] = decompToeplitz(T_row, 0, -1); % Get a standard against which to measure the accuracy of our decompToeplitz method. [TiInvT0, ldT0] = decompToeplitz(T_row, 0, 0); % Compute the inverse and log determinant using our decompToeplitz method. new_vals_times(3, 1) = norm(TiInvT(:) - TiInvT0(:), 2); % Compares accuracy to built-in Matlab algorithm. TiInvT = 0; % Clears some bulky variables for performance. clearvars is not used due to Octave incompatbility. TiInvT0 = 0; % Clears some bulky variables for performance. clearvars is not used due to Octave incompatbility. % Check the accuracy of the log determinant in our decompToeplitz method. new_vals_times(2, 1) = norm(ldT(:) - ldT0(:), 2); % Compares accuracy to built-in Matlab algorithm. % Check the accuracy of our traceToeplitz method. tr = traceToeplitz(A, T_row, T_row, -1); % Get a standard against which to measure the accuracy of our traceToeplitz method. trC = traceToeplitz(A, T_row, T_row, 1); % Compute the trace using our method. new_vals_times(4, 1) = norm(tr - trC, 2); % Compares accuracy to built-in Matlab algorithm. end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % test_iterToeplitz tests the number of iterations required for the preconditioned conjugate gradient method to converge using various preconditioners, and outputs the results in a column vector, new_vals_times. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [new_vals_times] = test_iterToeplitz(T_row, b) new_vals_times = zeros(3, 1); % Intialize new_vals_times as a vector of zeros. % Call solveToeplitzPCG and specify the unconditioned conjugate gradient method. [xPCG, k] = solveToeplitzPCG(T_row, b, 'None'); new_vals_times(1, 1) = k; % Stores the number of iterations required to converge. % Call solveToeplitzPCG and specify the Chan preconditioned conjugate gradient method. [xPCG, k] = solveToeplitzPCG(T_row, b, 'Chan'); new_vals_times(2, 1) = k; % Stores the number of iterations required to converge. % Call solveToeplitzPCG and specify the Circulant Embedding preconditioned conjugate gradient method. t = cputime; [xPCG, k] = solveToeplitzPCG(T_row, b, 'Circulant_Embedding'); new_vals_times(3, 1) = k; % Stores the number of iterations required to converge. end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % test_complexityToeplitz_small tests the complexity of built-in MATLAB and Toeblitz methods limited by O(n^2) storage, and outputs the results in a column vector, new_vals_times. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[new_vals_times] = test_complexityToeplitz_small(T_row, A, b); new_vals_times = zeros(8, 1); % Intialize new_vals_times as a vector of zeros. T = toeplitz(T_row); % Create the full matrix representation for built-in MATLAB methods. % Check the speed of our preconditioned conjugate gradient method. t = cputime; solveToeplitz(T_row, b); % Get the result using our method. new_vals_times(1, 1) = cputime - t; % Saves the time to our data matrix. % Check the speed of our logdetToeplitz method. t = cputime; logdetToeplitz(T_row); % Get the result using our method. new_vals_times(2, 1) = cputime - t; % Saves the time to our data matrix. % Check the speed of our invToeplitz method. t = cputime; invToeplitz(T_row); % Get the result using our method. new_vals_times(3, 1) = cputime - t; % Saves the time to our data matrix. % Check the speed of our traceToeplitz method. t = cputime; traceToeplitz(A, T_row, T_row, 1); % Compute the trace using our method. new_vals_times(4, 1) = cputime - t; % Saves the time to our data matrix. % Check the speed of the built-in matrix equation solver. t = cputime; mldivide(T, b); % Get the result using the built-in MATLAB method. new_vals_times(5, 1) = cputime - t; % Saves the time to our data matrix. % Check the speed of the built-in log determinant method. t = cputime; logdet(T); % Get the result using the built-in MATLAB method. new_vals_times(6, 1) = cputime - t; % Saves the time to our data matrix. % Check the speed of the built-in matrix inversion method. t = cputime; inv(T); % Get the result using the built-in MATLAB method. new_vals_times(7, 1) = cputime - t; % Saves the time to our data matrix. % Check the speed of the built-in matrix multiply and trace method. t = cputime; trace(A*T); % Get the result using the built-in MATLAB method. % sum(sum(A.*transpose(T))); % Get the result using the built-in MATLAB method. This line is actually a vast improvement over simply taking trace(A*T). new_vals_times(8, 1) = cputime - t; % Saves the time to our data matrix. end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % test_complexityToeplitz_large tests the complexity of Toeblitz methods that use O(n) storage, and outputs the results in a column vector, new_vals_times. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[new_vals_times] = test_complexityToeplitz_large(T_row, b); new_vals_times = zeros(2, 1); % Intialize new_vals_times as a vector of zeros. % Check the speed of our preconditioned conjugate gradient method. t = cputime; solveToeplitz(T_row, b); % Get the result using our method. new_vals_times(1, 1) = cputime - t; % Saves the time to our data matrix. % Check the speed of our logdetToeplitz method. t = cputime; logdetToeplitz(T_row); % Compute just the log determinant using our logdetToeplitz method. new_vals_times(2, 1) = cputime - t; % Saves the time to our data matrix. end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % grand_loglogplot plots the computational time used for calculations vs. matrix size for several methods on a log-log axis, allowing us to compare the complexities of the various methods. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [a_graph] = grand_loglogplot(sizes_vec_small, sizes_vec_large, data_mat_small, data_mat_large, data_labels, yaxis_label, save_flag, folder, file_name, runMode, fig_save) sizes_vec_large = horzcat(sizes_vec_small, sizes_vec_large); % Adds sizes_vec_small to the beginning of sizes_vec_large so that the full matrix size series in present. a_graph = figure; % Create a figure. loglog(sizes_vec_large, data_mat_large(1, 1:end), 'r', sizes_vec_large, data_mat_large(2, 1:end), 'b', sizes_vec_small, data_mat_small(1, 1:end), 'm', sizes_vec_small, data_mat_small(2, 1:end), 'g', sizes_vec_small, data_mat_small(3, 1:end), 'r--', sizes_vec_small, data_mat_small(4, 1:end), 'b--', sizes_vec_small, data_mat_small(5, 1:end), 'm--', sizes_vec_small, data_mat_small(6, 1:end), 'g--', 'LineWidth',3); % Plot our data. xlabel('Matrix Size', 'FontSize', 18); % Label the x-axis. ylabel(yaxis_label,'FontSize', 18); % Label the y-axis. set(gca,'fontsize', 18); % Set the general font size. legend(char(data_labels(1, 1)), char(data_labels(2, 1)), char(data_labels(3, 1)), char(data_labels(4, 1)), char(data_labels(5, 1)), char(data_labels(6, 1)), char(data_labels(7, 1)), char(data_labels(8, 1))) % Label the legend using the provided strings. % Save our plot if save_flag is on. Otherwise, display it. if save_flag print(a_graph, strcat(folder, file_name), '-dpdf'); if fig_save print(a_graph, strcat(folder, file_name), '-dfig'); end close all else a_graph; end %{ % View basic semilogy and linear plots of our data when in runMode == -1. This helps us see what's actually being plotted, since loglog has trouble with values at or near 0. Use for debugging. if runMode == -1 figure; semilogy(sizes_vec_large, data_mat_large(1, 1:end), 'r', sizes_vec_large, data_mat_large(2, 1:end), 'b', sizes_vec_small, data_mat_small(1, 1:end), 'm', sizes_vec_small, data_mat_small(2, 1:end), 'g', sizes_vec_small, data_mat_small(3, 1:end), 'r--', sizes_vec_small, data_mat_small(4, 1:end), 'b--', sizes_vec_small, data_mat_small(5, 1:end), 'm--', sizes_vec_small, data_mat_small(6, 1:end), 'g--', 'LineWidth',3); figure; plot(sizes_vec_large, data_mat_large(1, 1:end), 'r', sizes_vec_large, data_mat_large(2, 1:end), 'b', sizes_vec_small, data_mat_small(1, 1:end), 'm', sizes_vec_small, data_mat_small(2, 1:end), 'g', sizes_vec_small, data_mat_small(3, 1:end), 'r--', sizes_vec_small, data_mat_small(4, 1:end), 'b--', sizes_vec_small, data_mat_small(5, 1:end), 'm--', sizes_vec_small, data_mat_small(6, 1:end), 'g--', 'LineWidth',3); end %} end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % grand_regressionplot plots the log of the computational time used for calculations vs. the log of the matrix size for several methods on a linear axis, effectively producing the same figure as grand_loglogplot. However, this format allows us to use MATLAB's simple regression tools to determine the best fit slopes for each data series obtained. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [a_graph] = grand_regressionplot(sizes_vec_small, sizes_vec_large, data_mat_small, data_mat_large, data_labels, yaxis_label, save_flag, folder, file_name, runMode, fig_save) sizes_vec_large = horzcat(sizes_vec_small, sizes_vec_large); % Adds sizes_vec_small to the beginning of sizes_vec_large so that the full matrix size series in present. % Get the log of all the data, so that a normal plot will look the same as a loglog plot. sizes_vec_large = log(sizes_vec_large); sizes_vec_small = log(sizes_vec_small); data_mat_large = log(data_mat_large); data_mat_small = log(data_mat_small); a_graph = figure; % Create a figure. plot(sizes_vec_large, data_mat_large(1, 1:end), 'r', sizes_vec_large, data_mat_large(2, 1:end), 'b', sizes_vec_small, data_mat_small(1, 1:end), 'm', sizes_vec_small, data_mat_small(2, 1:end), 'g', sizes_vec_small, data_mat_small(3, 1:end), 'r--', sizes_vec_small, data_mat_small(4, 1:end), 'b--', sizes_vec_small, data_mat_small(5, 1:end), 'm--', sizes_vec_small, data_mat_small(6, 1:end), 'g--', 'LineWidth',3); % Plot our data. xlabel('Matrix Size', 'FontSize', 18); % Label the x-axis. ylabel(yaxis_label,'FontSize', 18); % Label the y-axis. set(gca,'fontsize', 18); % Set the general font size. legend(char(data_labels(1, 1)), char(data_labels(2, 1)), char(data_labels(3, 1)), char(data_labels(4, 1)), char(data_labels(5, 1)), char(data_labels(6, 1)), char(data_labels(7, 1)), char(data_labels(8, 1))) % Label the legend using the provided strings. % Save our plot if save_flag is on. Otherwise, display it. if save_flag print(a_graph, strcat(folder, file_name), '-dpdf'); if fig_save print(a_graph, strcat(folder, file_name), '-dfig'); end close all else a_graph; end %{ % View basic semilogy and linear plots of our data when in runMode. This helps us see what's actually being plotted, since loglog has trouble with values at or near 0. if runMode == -1 figure; semilogy(sizes_vec_large, data_mat_large(1, 1:end), 'r--', sizes_vec_large, data_mat_large(2, 1:end), 'g--', sizes_vec_small, data_mat_small(1, 1:end), 'b--', sizes_vec_small, data_mat_small(2, 1:end), 'y--', sizes_vec_small, data_mat_small(3, 1:end), 'r', sizes_vec_small, data_mat_small(4, 1:end), 'g', sizes_vec_small, data_mat_small(5, 1:end), 'b', sizes_vec_small, data_mat_small(6, 1:end), 'y', 'LineWidth',3); figure; plot(sizes_vec_large, data_mat_large(1, 1:end), 'r--', sizes_vec_large, data_mat_large(2, 1:end), 'g--', sizes_vec_small, data_mat_small(1, 1:end), 'b--', sizes_vec_small, data_mat_small(2, 1:end), 'y--', sizes_vec_small, data_mat_small(3, 1:end), 'r', sizes_vec_small, data_mat_small(4, 1:end), 'g', sizes_vec_small, data_mat_small(5, 1:end), 'b', sizes_vec_small, data_mat_small(6, 1:end), 'y', 'LineWidth',3); end %} end toeblitz-1.0/util/multiplyCirculant.m0000775000442500044270000000207512202534277021150 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This uses a fast Fourier transform to compute the product y = Cx, where C is a circulant matrix and x is a column vector of appropriate size. This is a simplified version of the Toeplitz matrix multiplication on page 202 of "Matrix Computations", 3rd edition by Gene H. Golub. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [y] = multiplyCirculant(C, x) n = size(C, 2); %Gets the dimension of the circulant matrix. embed_circulant = C(1:end, 1); %Extracts information from first column of circulant matrix. y = ifft(fft(embed_circulant).*fft(x)); %Computes the result of the circulant matrix multiplication in a fast way. endtoeblitz-1.0/util/convertToeplitz.m0000775000442500044270000000544512202534277020643 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This function automatically converts a Toeplitz matrix into a row vector containing all necessary information, or expands that row vector into a full Toeplitz matrix, depending on the input. For non-symmetric matrices, this function will output both the first row and first column. % % The calling syntax is % % T_row = convertToeplitz(T), where T is a symmetric Toeplitz matrix and T_row is its first row. % [T_row, T_col] = convertToeplitz(T), where T is an asymmetric Toeplitz matrix, T_row is its first row, and T_col is its first column. % T = convertToeplitz(T_row), where T_row is the first row of a symmetric Toeplitz matrix, and T is the full matrix. % T = convertToeplitz(T_row, T_col), where T_row is the first row of an asymmetric Toeplitz matrix T, T_col is the first column of that same matrix, and T is the full matrix. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [alt_form, alt_form2] = convertToeplitz(one_form, one_form2) if nargin == 1 % Fill in a value for one_form2 when it is not specified. one_form2 = 0; % Create a placeholder value for one_form2 when it is not specified in the input. This value is not used. end if size(one_form, 1) == size(one_form, 2) % This is the case in which the first input was a Toeplitz matrix. We only check that it is a square matrix, for the sake of efficiency. We leave it to the caller to make sure that the input is Toeplitz. alt_form = one_form(1, 1:end); if not(isequal(transpose(one_form(1:end, 1)), one_form(1, 1:end))) % This is the case in which the input was a non-symmetric Toeplitz matrix. alt_form2 = one_form(1:end, 1); % In this case we also output the first column of the original matrix as a second row in our output. end elseif size(one_form, 1) == 1 % This is the case in which the first input was a row vector. if nargin == 2 % If there were two input arguments, the Toeplitz matrix is not necessary symmetric. alt_form = toeplitz(one_form2, one_form); % We use the convention opposite that of the built-in toeplitz() method. So we reverse the order of our row and column vector to generate the correct Toeplitz matrix. elseif nargin == 1 % If there is only one input, generate a symmetric Toeplitz matrix. alt_form = toeplitz(one_form); % Generate a symmetric Toeplitz matrix. else error('Too many input arguments!') end else error('Invalid input!') end endtoeblitz-1.0/util/startup_checkMEX.m0000775000442500044270000000474312202534277020641 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang and John Cunningham, 2013 % The following code checks for the relevant MEX files (such as .mexa64 % or .mexglx, depending on the machine architecture), and it creates the % mex file if it can not find the right one. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [] = startup_checkMEX(file_directory, file_name) toeplitz_directory = sprintf('%s', pwd); % Finds the Toeplitz directory in a stable way. toeplitz_directory = strrep(toeplitz_directory, '\', '/'); % Replaces Windows backslashes with forward slashes. file_directory = [toeplitz_directory '/' file_directory]; % Creates the full file directory. if not(exist(sprintf([file_directory '/' file_name '.' mexext]), 'file')) % Check if the relevant files exist. try % If not, try to create them. if isOctave eval(sprintf(['mkoctfile --mex --output ' file_directory '/' file_name '.mex ' file_directory '/' file_name '.c'])) else eval(sprintf(['mex -outdir ' file_directory ' ' file_directory '/' file_name '.c'])); % This line actually compiles the C code and creates the MEX file. end fprintf(['NOTE: the relevant ' file_name ' mex files were not found. They have been created.\n']); % Success message. catch % If the files cannot be created, display the following error message. fprintf(['NOTE: the relevant ' file_name ' mex files were not found, and your machine failed to create them.\n']); fprintf(' This usually means that you do not have the proper C/MEX compiler setup.\n'); fprintf(' The code will still run identically, albeit slower (perhaps considerably).\n'); fprintf(' Please read the README file section on the use of C/MEX.\n'); end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The function isOctave detects whether the user is using Octave, returning true if this is the case, and false if it is not (meaning the user is using MATLAB, most likely). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [OctaveMode] = isOctave() OctaveMode = exist('OCTAVE_VERSION') ~= 0; end toeblitz-1.0/testToeplitz.m0000775000442500044270000002664112202535720017160 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, John P Cunningham % This tests a toeplitz function with different sizes of matrices and times it to evaluate performance and accuracy. It takes two required input arguments, method and sizes_vec, and two optional arguments, runs and save_flag. The calling syntax is % testToeplitz(method, sizes_vec[, runs, save_flag]). % If three arguments are given, the code will assume that the third argument is runs. To specify save_flag, runs must also be specified. % % Valid inputs: % method = 'decompToeplitz', 'traceToeplitz', 'solveToeplitz', or 'solveToeplitzPCG' % sizes = a row vector of positive integer sizes over which to evaluate the performance of that function; for example, [5:10] will cause testToeplitz to produce results for matrices at sizes 5, 6, 7, 8, 9, and 10, and display relevant plots. % runs = a positive integer which indicates the number of runs over which to average results for each matrix size and algorithm. More runs will give smoother results and nicer graphs, but will be more computationally costly. This defaults to 3 if not specified. % save_flag = a boolean value, 0 or 1, which indicates whether plots should be displayed (0) or saved to the ./results directory. This defaults to 0 (only display, no saving) if not specified. % % Some example calls: % testToeplitz('decompToeplitz', [2:10:202]) % testToeplitz('traceToeplitz', [500:600]) % testToeplitz('solveToeplitz', [5 10 100 10^3 10^4], 10) % testToeplitz('solveToeplitzPCG', [400:100:1000], 1, 1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This is the main function. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [] = testToeplitz(method, sizes_vec, runs, save_flag) % The following code sets default values for runs and save_flag if they are not specified. if nargin < 4 % Check to see if any arguments are missing. save_flag = 0; % Default value for save_flag is false. if nargin < 3 runs = 3; % Default number of runs to average over. end end % The following code finds a folder to save the result plots in, if save_flag is set to 1. folder = '0'; if save_flag for i = 1:999 % Searches from 001 to 999 for an unused folder name of the form 'test_XXX', where XXX is the three-digit number. if folder == '0' % Runs the following code if we have not yet found a suitable folder. test_folder = sprintf('./results/test_%03d/', i); % Finds the subdirectory in which the files will be placed. if not(exist(test_folder)) % Checks to see if there is already a folder with that name in the ./results/ directory. If not, it creates the folder and proceeds. folder = test_folder; % Changes the value of folder, stopping future checks for folder names. mkdir(folder); % Actually creates the directory. fprintf(['Printing to directory ' folder '.\n']) % Informs the user where the output will be saved. end end end end % The following code tests the accuracy and speed of our traceToeplitz method. Routines for testing other functions work similarly, and are less extensively commented. if strcmp(method, 'traceToeplitz') % This checks if 'traceToeplitz' is the method specified for testing. trace_strs = cell(4, 1); % This creates a cell array of the names of our different methods for use in labeling plots and legends. trace_strs(1, 1) = cellstr('trace(A*T)'); trace_strs(2, 1) = cellstr('traceToeplitz (Automated Choice)'); trace_strs(3, 1) = cellstr('traceToeplitz (C/MEX)'); trace_strs(4, 1) = cellstr('traceToeplitz (Matlab)'); vals_times = zeros(8, length(sizes_vec)); % This initializes the matrix in which all our test results will be (temporarily) stored. for i = 1:length(sizes_vec) % This loops over each the sizes of matrices we want to test. current_size_str = num2str(sizes_vec(i)); % Finds the current matrix size being worked on. fprintf(['Working on matrix size ' current_size_str '.\n']) % Informs the user what matrix size we're up to. new_vals_times = zeros(8, 1); % This initializes new_vals_times (and resets it for each size of matrix tested). for j = 1:runs % For each size matrix, this loops 'runs' times to obtain an average result. run_str = num2str(j); % Finds the current run being worked on. fprintf(['Now on run number ' run_str '.\n']) % Informs the user which run we are on. T_row = genPDToeplitz(sizes_vec(i)); % This produces the first row of a (more or less) random positive definite Toeplitz matrix. new_vals_times = new_vals_times + test_traceToeplitz(T_row); % This calculates the trace of the product with an arbitrary matrix in various ways, and adds the timings and accuracy metrics to a cumulative total in new_vals_times (see test_traceToeplitz documentation for more detail). end vals_times(1:8, i) = new_vals_times(1:8, 1)/runs; % This saves the results in new_vals_times to the appropriate column in vals_times, dividing by the number of runs to obtain an average. end % Show or save results % These are results for computational time. test_plot('all', sizes_vec, vals_times(2:2:8, 1:end), trace_strs, 'Computation time (s)', save_flag, folder, 'test_trace_time', 0); % These are results for accuracy. test_plot('linear', sizes_vec, vals_times(1:2:7, 1:end), trace_strs, 'Norm difference vs. trace(A*T)', save_flag, folder, 'test_trace_norm', 0); elseif strcmp(method, 'decompToeplitz') decomp_strs = cell(4, 1); decomp_strs(1, 1) = cellstr('decompToeplitz (Automated Choice)'); decomp_strs(2, 1) = cellstr('inv(T)'); decomp_strs(3, 1) = cellstr('C/MEX Zohar'); decomp_strs(4, 1) = cellstr('Matlab Zohar'); vals_times = zeros(12, length(sizes_vec)); for i = 1:length(sizes_vec) current_size_str = num2str(sizes_vec(i)); % Finds the current matrix size being worked on. fprintf(['Working on matrix size ' current_size_str '.\n']) % Informs the user what matrix size we're up to. new_vals_times = zeros(12, 1); for j = 1:runs run_str = num2str(j); % Finds the current run being worked on. fprintf(['Now on run number ' run_str '.\n']) % Informs the user which run we are on. T_row = genPDToeplitz(sizes_vec(i)); % Make random positive definite symmetric toeplitz matrix. % Mow calculate the inverse in various ways, and save the timings and accuracy metrics (see test_decompToeplitz documentation for more detail). new_vals_times = new_vals_times + test_decompToeplitz(T_row, i); end vals_times(1:12, i) = new_vals_times(1:12, 1)/runs; % Update our array of results with the average of all trials, after each size matrix in sizes_vec is tested. end % Show or save results % These are results for computational time. test_plot('all', sizes_vec, vals_times(3:3:12, 1:end), decomp_strs, 'Computation time (s)', save_flag, folder, 'test_decomp_time', 0); % These are results for accuracy. test_plot('linear', sizes_vec, vals_times(1:3:10, 1:end), decomp_strs, 'Frobenius norm difference vs. inv(T)', save_flag, folder, 'test_decomp_invnorm', 0); test_plot('linear', sizes_vec, vals_times(2:3:11, 1:end), decomp_strs, 'Norm difference vs. logdet(T)', save_flag, folder, 'test_decomp_logdetnorm', 0); elseif strcmp(method, 'solveToeplitz') solve_strs = cell(4, 1); solve_strs(1, 1) = cellstr('Matlab mldivide'); solve_strs(2, 1) = cellstr('Levinson-Durbin'); solve_strs(3, 1) = cellstr('Unconditioned Conjugate Gradient'); solve_strs(4, 1) = cellstr('Preconditioned Conjugate Gradient'); vals_times = zeros(8, length(sizes_vec)); for i = 1:length(sizes_vec) current_size_str = num2str(sizes_vec(i)); % Finds the current matrix size being worked on. fprintf(['Working on matrix size ' current_size_str '.\n']) % Informs the user what matrix size we're up to. new_vals_times = zeros(8, 1); for j = 1:runs run_str = num2str(j); % Finds the current run being worked on. fprintf(['Now on run number ' run_str '.\n']) % Informs the user which run we are on. % Make a random positive definite symmetric toeplitz matrix. T_row = genPDToeplitz(sizes_vec(i)); % Generates an arbitrary vector b to use in the equation Tx = b when solving for x. b = [1:sizes_vec(i)]; b = b(:); % Now solve the system Tx = b in various ways, and return the timings and accuracy metrics. new_vals_times = new_vals_times + test_solveToeplitz(T_row, b); %Tests the accuracy and speed of our algorithms for this size matrix. end vals_times(1:8, i) = new_vals_times(1:8, 1)/runs; % Adds the results from new_vals_times to val_times. end % Show or save results %These are results for computational time. test_plot('all', sizes_vec, vals_times(2:2:8, 1:end), solve_strs, 'Computation time (s)', save_flag, folder, 'test_solve_time', 0); %These are results for accuracy. test_plot('linear', sizes_vec, vals_times(1:2:7, 1:end), solve_strs, 'Norm difference vs. Built-in mldivide', save_flag, folder, 'test_solve_norm', 0); elseif strcmp(method, 'solveToeplitzPCG') PCG_strs = cell(4, 1); PCG_strs(1, 1) = cellstr('Unconditioned Conjugate Gradient'); PCG_strs(2, 1) = cellstr('MATLAB pcg'); PCG_strs(3, 1) = cellstr('Chan Preconditioner'); PCG_strs(4, 1) = cellstr('Circulant Embedding Preconditioner'); vals_times = zeros(12, length(sizes_vec)); for i = 1:length(sizes_vec) current_size_str = num2str(sizes_vec(i)); % Finds the current matrix size being worked on. fprintf(['Working on matrix size ' current_size_str '.\n']) % Informs the user what matrix size we're up to. new_vals_times = zeros(12, 1); for j = 1:runs run_str = num2str(j); % Finds the current run being worked on. fprintf(['Now on run number ' run_str '.\n']) % Informs the user which run we are on. % Make random positive definite symmetric toeplitz matrix. T_row = genPDToeplitz(sizes_vec(i)); % Generates an arbitrary vector b to use in the equation Tx = b when solving for x. b = [1:sizes_vec(i)]; b = b(:); % Now solve the system Tx = b in various ways, and return the timings and accuracy metrics. new_vals_times = new_vals_times + test_solveToeplitzPCG(T_row, b); %Tests the accuracy and speed of our algorithms for this size matrix. end vals_times(1:12, i) = new_vals_times(1:12, 1)/runs; % Adds the results from new_vals_times to val_times. end % Show or save results % These are results for computational time. test_plot('all', sizes_vec, vals_times(2:3:11, 1:end), PCG_strs, 'Computation time (s)', save_flag, folder, 'test_PCG_time', 0); % These are results for accuracy. test_plot('linear', sizes_vec, vals_times(1:3:10, 1:end), PCG_strs, 'Norm difference vs. Built-in mldivide', save_flag, folder, 'test_PCG_norm', 0); %These are results for the number of iterations. test_plot('linear', sizes_vec, vals_times(3:3:12, 1:end), PCG_strs, 'Number of Iterations', save_flag, folder, 'test_PCG_iterations', 0); else error('No such method found.') end end toeblitz-1.0/invToeplitz.m0000775000442500044270000000275512202534277017003 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This is simply a wrapper function that calls decompToeplitz in a way that produces both ld, the log determinant of T, and Ti, the inverse. Note that the log determinant is computed essentially for free when the inverse is also obtained. We can specify a runMode here to run decompToeplitz in the desired mode. Here are the accepted values of runMode: % runMode = 0 % This calls the C/MEX implementation of Zohar's modification of the Trench algorithm. % runMode = 1 % This calls the backup Matlab implementation of Zohar's modification of the Trench algorithm. % runMode = -1 % This calls the generic MATLAB function logdet. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[Ti, ld] = invToeplitz(T_row, runMode) if nargin < 2 % Handles the case when runMode is not specified. [ld, Ti] = decompToeplitz(T_row, 0); % This allows decompToeplitz to automatically choose the best method. else % Handles the case when runMode is specified. [ld, Ti] = decompToeplitz(T_row, 0, runMode); % This forces decompToeplitz to use the specified method. end end toeblitz-1.0/decomp/0000775000442500044270000000000012202534277015531 5ustar cunninghamcunninghamtoeblitz-1.0/decomp/logdetToeplitzFastZohar.mex0000775000442500044270000002163612202534277023112 0ustar cunninghamcunninghamELF44 ($!T T HH$$QtdRtdttGNU\y H|20Ax[`'nt     ` |CEqX4K + |"tX_0 {( ( ?A 0 H  __gmon_start___init_fini__cxa_finalize_Jv_RegisterClassesmexFunctionmexErrMsgTxtmxGetNmxCreateDoubleScalarmxGetPrmxMalloclogmxFreeliboctinterp.soliboctave.solibcruft.soliblapack.so.3gflibblas.so.3gflibfftw3.so.3libfftw3f.so.3libreadline.so.6libncurses.so.5libdl.so.2libhdf5.so.6libz.so.1libgfortran.so.3libstdc++.so.6libm.so.6libgcc_s.so.1libc.so.6_edata__bss_start_endGLIBC_2.0GLIBC_2.1.3Y ii qsi $           US[øt.)X[ hhhhh h($h0(h8p,h@`UVS4u]t0$l8)9s 889rƃ4[^]US.ktt $Ѓ[]Ë$ÐUWVSlu#}tp$$E$EU $E$}<$EYEs<$FEs<$3t<$!Et<$EUM11}ĉu܍v4׉19EwEċUԋMuȃMثX}1E}ȍv}ЋEЋM}܋}}E1t& ƃ9~؋M19~E܃EM0E1UЉUU2U \;EuȋM܋}̃E9EUMLE}t5؉u܋u11Uԃ$]M9Ew܋u܋}$]E1҉UUEEmE$M $4$}ԉ<$Eĉ$l[^_]Ã}$}UċMEثXPE4$$z$z$y $tyUVSt&Ћu[^]US[àY[1 input arg required.2 output args required.Inadequate space on heap for r.Inadequate space on heap for lambda.Inadequate space on heap for ghat.Inadequate space on heap for ghatold.Inadequate space on heap for gamma.?"/9JYcq 0 H op  H ooofov$ GCC: (Ubuntu 4.4.3-4ubuntu5.1) 4.4.3AmexFunction}  int25lgbV1 ]1(,:1]X 1(D1iP4z  4zu r4z I4zuT ?4z 4z y4zu n5D i6(/ j6(n c = o u= 6% $ > $ > : ; $ > .? : ; ' @: ; I4: ; I 4: ; I 4: ; I 4: ; I 4: ; I  I&I! /home/cunni/toeplitz/wemmick/toeplitz/decomp/usr/include/octave-3.2.3/octavelogdetToeplitzFastZohar.cmex.h1n *@8@@Y-=k-=k-/k-=k-=m>,5jbNbic&(dZH dz?9[9[^dfd*f;Y-/:>Lr^x8$;.gs,Gh.444| AAB Inlhslogdetnrhsunsigned charshort unsigned intmxArrayplhsghatprhslambdaT_rowmexFunctionlong long unsigned intgammaghatoldlong long intGNU C 4.4.3/home/cunni/toeplitz/wemmick/toeplitz/decomp/logdetToeplitzFastZohar.cshort intlong doublettAuuAu Au !!Au! !VuVuVAuQ/VP+RR2W2?R\P.symtab.strtab.shstrtab.note.gnu.build-id.gnu.hash.dynsym.dynstr.gnu.version.gnu.version_r.rel.dyn.rel.plt.init.text.fini.rodata.eh_frame.ctors.dtors.jcr.dynamic.got.got.plt.data.bss.comment.debug_aranges.debug_pubnames.debug_info.debug_abbrev.debug_line.debug_frame.debug_str.debug_loc$2X.opp<8 @Hoff"Uo@d  m H v000q``|HH H d d P P H 0$ $( (0(%M m"%4*0$5"@@#4 !pf 0 `  H  d P $ (   (5 K( Z, h tP    $ H  !.( :0 ? N( U ir"A 0 crtstuff.c__CTOR_LIST____DTOR_LIST____JCR_LIST____do_global_dtors_auxcompleted.7021dtor_idx.7023frame_dummy__CTOR_END____FRAME_END____JCR_END____do_global_ctors_auxlogdetToeplitzFastZohar.c__i686.get_pc_thunk.bx_GLOBAL_OFFSET_TABLE__DYNAMIC__DTOR_END____dso_handle_finimexErrMsgTxt__bss_start_end__gmon_start___edata_Jv_RegisterClassesmxMalloc__cxa_finalize@@GLIBC_2.1.3mexFunctionmxGetPr_initlog@@GLIBC_2.0mxGetNmxCreateDoubleScalarmxFreetoeblitz-1.0/decomp/logdet.m0000775000442500044270000000071212202534277017170 0ustar cunninghamcunningham% log(det(A)) where A is positive-definite. % This is faster and more stable than using log(det(A)). % Written by Tom Minka % (c) Microsoft Corporation. All rights reserved. % Additional comments added by William Zhang, 2013. function [y] = logdet(A) U = chol(A); % Standard Cholesky decomposition A = U'*U y = 2*sum(log(diag(U))); % The product of the diagonal % of any triangular matrix is the determinant. endtoeblitz-1.0/decomp/decompToeplitzFastZohar.c0000775000442500044270000001523612202534277022533 0ustar cunninghamcunningham/*================================================================= * Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices * Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) *=================================================================*/ /*================================================================= * * decompToeplitzFastZohar.c completes the inversion of a symmetric * positive definite Toeplitz matrix. This * function is a subroutine of the decompToeplitz.m, which follows the * algorithm of Zohar 1969, a modification of the algorithm of W.F. Trench. * * The calling syntax is: * * [ Ti , logdetT ] = decompToeplitzFastZohar(T_row) * * T_row is the first row of the positive definite symmetric Toeplitz matrix T of dimension n+1 * (following the convention of p599 in the Zohar paper) * Ti is the inverse of T, same as inv(T) in MATLAB. * logdet is the log determinant of T, NOT of Ti. * * NOTE: This code is used to speed up the Toeplitz inversion as much * as possible. Accordingly, no error checking is done. The onus is * on the caller (which should be decompToeplitz.m) to pass the correct arguments. * * NOTE: This is superior to the Trench algorithm in Golub and * Van Loan's "Matrix Computations", since the Zohar version * also computes the log determinant for free, which is essential in the * application for which this algorithm was coded. * * John P Cunningham * 2009 * Modified by William Zhang, 2013 *=================================================================*/ /* Required Packages */ #include #include /* Input Arguments */ #define T_row_IN prhs[0] /* Output Arguments */ #define Ti_OUT plhs[0] #define LD_OUT plhs[1] void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] ) { /* Declare variables */ double *T_row, *Ti, *logdet, *r, *lambda, *ghat, *ghatold, *gamma; unsigned int n; int i,j; /* Check for proper number of arguments*/ if (nrhs != 1) { mexErrMsgTxt("1 input arg required."); } else if (nlhs > 2) { mexErrMsgTxt("2 output args required."); } /* Check the dimensions of the arguments */ /* This expects T to be an n+1 by n+1 matrix */ n = mxGetN(T_row_IN)-1; /* Create a matrix for the return argument */ Ti_OUT = mxCreateDoubleMatrix(n+1, n+1, mxREAL); LD_OUT = mxCreateDoubleScalar(0); /* LD_OUT is a pointer to a 1x1 mxArray double initialized to 0. */ /* Assign pointers to some input and output parameters */ Ti = mxGetPr(Ti_OUT); logdet = mxGetPr(LD_OUT); T_row = mxGetPr(T_row_IN); /*Allocate space on heap for various variables.*/ r = (double*) mxMalloc(n*sizeof(double)); if (r == NULL) { mexErrMsgTxt("Inadequate space on heap for r."); } lambda = (double*) mxMalloc(n*sizeof(double)); if (lambda == NULL) { mexErrMsgTxt("Inadequate space on heap for lambda."); } ghat = (double*) mxMalloc(n*sizeof(double)); if (ghat == NULL) { mexErrMsgTxt("Inadequate space on heap for ghat."); } ghatold = (double*) mxMalloc(n*sizeof(double)); if (ghatold == NULL) { mexErrMsgTxt("Inadequate space on heap for ghatold."); } gamma = (double*) mxMalloc(n*sizeof(double)); if (gamma == NULL) { mexErrMsgTxt("Inadequate space on heap for gamma."); } /* Define r, the normalized row from T_row(1,2) to T_row(1,n+1) */ for (i = 0; i < n ; i++) { r[i] = T_row[0*(n+1) + (i+1)]/T_row[0*(n+1) + 0]; } /* Initialize the values of lambda and g for the algorithm */ lambda[0] = 1 - pow(r[0],2); ghat[0] = -r[0]; /* Recursion to build g and lambda */ for (i=0 ; i < n-1 ; i++) { /* Calculate gamma */ gamma[i] = -r[i+1]; /* Note that ghat, g_i, etc are i+1 elements long */ for (j=0 ; j < i+1 ; j++) { gamma[i] -= r[j]*ghat[j]; } /* Calculate next ghat */ /* First store ghatold, which is i+1 elts long */ for (j=0 ; j < i+1 ; j++) { ghatold[j] = ghat[j]; } /* The first element of the new ghat is derived from gamma and lambda. */ ghat[0] = gamma[i]/lambda[i]; /* The rest of the new ghat is formed from the old ghat, ghatold, as described on page 599 of the Zohar paper. */ for (j=1 ; j < i+2 ; j++) { ghat[j] = ghatold[j-1] + gamma[i]/lambda[i]*ghatold[i+1-j]; } /* Calculate lambda */ lambda[i+1] = lambda[i] - pow(gamma[i],2)/lambda[i]; } /* There are n lambdas, n g values, and n-1 gammas */ /* Note that n is not the dimension of the matrix, but one less. */ /* This corresponds with Zohar notation. */ /* Evaluate the matrix Ti (B_{n+1} in Zohar) */ /* NOTE ON MEX MATRIX INDEXING */ /* indexing is always (colidx * colLen + rowidx) */ /* Assign first upper left element */ Ti[ 0*(n+1) + 0 ] = 1/lambda[n-1]; /* Assign the first row and column*/ for (i = 1; i < n+1; i++) { /* First row */ Ti[ 0*(n+1) + i] = ghat[ n-1-(i-1) ]/lambda[ n-1 ]; /* First column */ Ti[ i*(n+1) + 0] = ghat[ n-1-(i-1) ]/lambda[ n-1 ]; } /* Fill in last row and column by persymmetry and symmetry with first row and column.*/ for (i = 1; i < n ; i++) { /* Last row */ Ti[ n*(n+1) + i ] = ghat[ n-1-(n-i-1) ]/lambda[ n-1 ]; /* Last column */ Ti[ i*(n+1) + n ] = ghat[ n-1-(n-i-1) ]/lambda[ n-1 ]; } /* Assign the last (bottom right) element by persymmetry */ Ti[ n*(n+1) + n ] = Ti[ 0*(n+1) + 0]; /* Fill in the interior of Ti_OUT */ for (i = 0; i < n/2 ; i++) { for (j = i; j < n-i-1; j++) { /* Calculate the value using equations on page 599 of Zohar. */ Ti[ (j+1)*(n+1) + (i+1) ] = (Ti[ j*(n+1) + i ] + (1/lambda[n-1])*(ghat[n-1-i]*ghat[n-1-j] - ghat[n-1-(n-1-i)]*ghat[n-1-(n-1-j)])); /* Use symmetry */ Ti[ (i+1)*(n+1) + (j+1) ] = Ti[ ((j+1)*(n+1)) + (i+1) ]; /* Use persymmetry: Recall there are n+1 elements, so 0<->n, 1<->n-1, and so on */ Ti[ (n-(i+1))*(n+1) + (n-(j+1)) ] = Ti[ ((j+1)*(n+1)) + (i+1) ]; Ti[ (n-(j+1))*(n+1) + (n-(i+1)) ] = Ti[ ((j+1)*(n+1)) + (i+1) ]; } } /* Normalize the entire matrix by T(1,1) so it is the properly scaled inverse */ for (i = 0; i < n+1 ; i++) { for (j = 0; j < n+1 ; j++) { Ti[ j*(n+1) + i ] = Ti[ j*(n+1) + i ]/T_row[ 0*(n+1) + 0]; } } /* Calculate the log determinant for free (essentially) by summing the logs of the lambdas. */ logdet[0] = 0; for (i=0 ; i < n; i++) { logdet[0] += log(lambda[i]); } /* Renormalize based on T_row(1,1) */ logdet[0] += (n+1)*log(T_row[ 0*(n+1) + 0 ]); /* Free allocated arrays */ mxFree(gamma); mxFree(ghatold); mxFree(ghat); mxFree(lambda); mxFree(r); return; }toeblitz-1.0/decomp/decompToeplitzFast.m0000775000442500044270000000564312202534277021542 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % John P Cunningham, 2009 % modified by William Zhang, 2013 % % decompToeplitzFast() % % This function is simply a wrapper for the C-MEX % function decompToeplitzFastZohar.mexa64 (or .mexglx, etc), % which is just a compiled version of invToeplitzFastZohar.c, % which should also be in this folder. Please see % that code for details. % % This algorithm inverts a positive definite (symmetric) % Toeplitz matrix in O(n^2) time, which is considerably better % than the O(n^3) that inv() offers. This follows the Trench % algorithm implementation of Zohar 1969. See that paper for % all explanation, as the C code is just an implementation of % the algorithm of p599 of that paper. % % This function also computes the log determinant, which is % calculated essentially for free as part of the calculation of % the inverse. This is often useful in applications when one really % needs to represent the inverse of a Toeplitz matrix. % % This function should be called from within decompToeplitz.m, % which adds a try block so that it can default to a native % MATLAB inversion (either inv() or a vectorized version of % the Zohar/Trench algorithm, depending on the size of the matrix) % should the MEX interface not work. This will happen, for example, % if you move to a new architecture and do not compile for .mexmaci % or similar (see mexext('all') and help mex for some info on this). % % Inputs: % T_row the first row of a positive definite (symmetric) Toeplitz matrix T, which % does NOT need to be scaled to be 1 on the main diagonal. % % Outputs: % Ti the inverse of T % ld the log determinant of T, NOT Ti. % % % NOTE: This code is used to speed up the Toeplitz inversion as much % as possible. Accordingly, no error checking is done. The onus is % on the caller (which should be decompToeplitz.m) to pass the correct args. % % NOTE: Whenever possible, do not actually invert a matrix. This code is % written just in case you really need to do so. Otherwise, for example % if you just want to solve inv(T)*x for some vector x, you are better off % using a fast inversion method, like PCG with fast matrix multiplication (used in solveToeplitz.m), % which could be something like an FFT method for the Toeplitz matrix. To % learn about this, see Cunningham, Sahani, Shenoy (2008), ICML, "Fast Gaussian % process methods for point process intensity estimation." %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Ti, ld] = decompToeplitzFast(T_row) [Ti, ld] = decompToeplitzFastZohar(T_row); end toeblitz-1.0/decomp/decompToeplitzFastZohar.mexw640000775000442500044270000002600012202534277023432 0ustar cunninghamcunninghamMZ@ !L!This program cannot be run in DOS mode. $`wwwwww܃wywvww݃wwwwwRichwPEd R"  #cJ@7]3d`P p 0.text `.rdata0@@.data`@$@.pdata P&@@.rsrc`(@@.reloc p*@BHSUVWATAWHLhLp)xD)@IHD)HAt H ! ~ H !<HM-E3DhL$Au֋Ή$XHfWHHGHOHHD$8HMH$EL$IHIHD$`HH$Hu H !IuLHu H IYLHD$XHu H ^I8H$Hu H @ =ILH$Hu H E3E3A|oAEHHMH+DD LAH I^ED I^ML A^ED I^ML uE;s*LAJ L+A+ADHH^EAuD 4fA(^ =AUf(HD$0A\ȉT$(A $DfAWA`H$H$I)$IHLH+I+MHD$@L+AHL$PL+HD$HfffHD9E3BT(3fAWAAAIHKH+DDHfDD H IYA\L YI\D YA\L YI\uE;},MAH L+A+BHHYA\uE3E3A|[AAIHNH+DDLHD H IHAHD HAHD HAHD HAuE;}!MAJ L+A+JHHHAuEQHDHD$ D$`A^>AD;AHD$PLIOHT8AAM+DDHHD$ H H IA^>YB0AXDAA^ >YJ(AXLIA^>YB AXDAA^ >YJAXLIuE;}GI HD$0LH+D$ EM+HTE+HHIA^>YBBXD AuA4>fA( HL$PHT$@^HD$0D$`HD9HHL$H\tH\$8$XHl$`($L$T$(AD($D($f(EMcIA^F!BA$`CD-E++FKTD$`GDmLS$XH ffBA(I H C^DABJ0C^L B(C0C^DABJ(C^L B A0C^DABJ C^LB BEHC^DABJC^L SHl$`L|$XD$XD;sMILˋH+AI NjA+DIHC^DA@IC^L IuIAAE*BADƍHDPCD-Eʼn$XA+AEE+IWD$XD HHD$ BA@H C^DJA*C^L BC2C^DBJC^L BA@C^DBA@AC^DBJC^L BA2EHC^DJHl$`L|$XHL$ E;sIDMAAEAC AIC^DA@C^DE;rHAME3E3LL$PAHA艄$DFKTID$XHT$0L|$ AAMA+LL$HɉL$(D; A+ƃL|$0DQƉD$|BD$lHEAHD$8FDT$tABT0IGHD$@AF$`+ȍAQƉ$֍BT$p$AƉ$ȉ$AD$xA@D$`D$hAA+I HD$8HL$HHL$@@QD$xDT$hDL$|AAYYEf(\DD+EC^LAIYf(XBEŋD+ABHHA@A+BL(HӋ$`HHD$@AIC^LHD$8@AYYE$\Yf(XBHHA@A+BL(HӋ$`HHD$@AIC^LPHD$8AYYEB\YXË$‹ABHHA@A+f(σBL(HӋ$`HHD$@AIC^LPHD$8@AYYE$\YXË$‹ABHHA@D$`DD$lA+BL(HӋT$pHHD$8HL$@H H D$`HD$8HL$@;|$tL|$XL$(HT$0D$XLL$PLT$H;l$(O Ɖ$`GDFT0IDEH+D$HM|AHD$ f(IIC^L$XYAYAHT$0\YBXDBË$`+BL(JHA+A+BL(JDH;nL|$XHl$ D$XLL$PDIHHAD$XLL$PHT$0Hl$ D;$Hl$`($E3ۅfD3CD-AABTD+ЍFE+F4D<B^EB ^EӍ2^EAI^Eu;s"΋A+^EHuAD;?H$L$EL$PHt*L$IHIXuL|$XEH$fH*YX9H$,I$IH$HA_A\_^][%%%%r%d%V%@%%@SH HHH H HuC#H#H H 3H [HHXHhHxL` AUAVAWH 3ML8#DeH%0HXH;t3HuAtcH LHH MLHHI;rZH9}tH9EtHMTHHEH V8H AH(L;uL;tLLIEGHH=EH=3eH%0HXH;t3Hut ;>H|H ekuHCH 4Eu HH?H9=Pt!H GtMĺI-H\$@Hl$HH|$PLd$XH A_A^A]HHXHpHxATH0ILXu9u 3ۉXtu7HHtЋ؉D$ LƋI0؉D$ LƋI؉D$ u5u1L3IL3IL9Mt L3IAӅtu7LƋI#ˋىL$ tHHtLƋIЋ؉D$ 3ۉ\$ H\$@Ht$HH|$PH0A\H\$Ht$WH IHuoLNjHH\$0Ht$8H _H M@SH HH `B HD$8Hu H ~H 2 HD$8H  HD$@H HLD$@HT$8HHL$8 HHL$@ HWHH [H(GHH(H\$WH H H= HHtHH;rH\$0H _H\$WH H H=x HHtHH;rH\$0H _HMZf9t3HcH6F6V6d6|66666606555555x5?@Inadequate space on heap for gamma.Inadequate space on heap for ghatold.Inadequate space on heap for ghat.Inadequate space on heap for lambda.Inadequate space on heap for r.2 output args required.1 input arg required.!3!2!2!2!h2! 2!x2!"*3# p`P020  t T 42t d 4R%"i#&i#}"o# 'd42 p20%$}$:'  4 2p2P B%%%`'%B  4 2p@550055046P0`470v7`7J7:7 7766666$65>6F6V6d6|66666606555555x5mxFree5mxMallocmxGetPrQmxCreateDoubleScalarOmxCreateDoubleMatrix_700mxGetNlibmx.dll1mexErrMsgTxtlibmex.dllpowlogMSVCR100.dll_malloc_crt_initterm_initterm_ecfree_encoded_null_amsg_exit__C_specific_handler__CppXcptFilter@__clean_type_info_names_internal[_unlockH__dllonexit_lock_onexitEncodePointerDecodePointerSleepDisableThreadLibraryCallsQueryPerformanceCounterGetTickCountGetCurrentThreadIdGetCurrentProcessIdGetSystemTimeAsFileTimeKERNEL32.dll R77777decompToeplitzFastZohar.mexw64mexFunctionu2-+] f32222222fx2$3U",3X"#H3##3#$3$$3$$3$%3%%3&#&3$&&3&'3 ':'3:'U'3`''30 HX`Z PAPADDINGXXPADDINGPADDINGXXPADDINGPADDINGXXPADDINGPADDINGXXPADDINGPADDINGXXPAD0 0toeblitz-1.0/decomp/logdetToeplitzFastZohar.c0000775000442500044270000001210612202534277022533 0ustar cunninghamcunningham/*================================================================= * Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices * Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) *=================================================================*/ /*================================================================= * * Adapted from decompToeplitzFastZohar.c. * logdetToeplitzFastZohar uses the Zohar algorithm to determine * the log determinant of a positive definite Toeplitz matrix, without * performing the computationally costly task of constructing the inverse. This * function is a subroutine of the logdetToeplitz.m, which follows the * algorithm of Zohar 1969, a modification of the algorithm of W.F. Trench. * * The calling syntax is: * * [ logdetT ] = logdetToeplitzFastZohar(T_row) * * T_row is the first row of the positive definite symmetric Toeplitz matrix T of dimension n+1 * (following the convention of p599 in the Zohar paper) * logdet is the log determinant of T. * * NOTE: This code is used to speed up the Toeplitz inversion as much * as possible. Accordingly, no error checking is done. The onus is * on the caller (which should be decompToeplitz.m) to pass the correct arguments. * * NOTE: This is superior to the Trench algorithm in Golub and * Van Loan's "Matrix Computations", since the Zohar version * also computes the log determinant for free, which is essential in the * application for which this algorithm was coded. * * John P Cunningham * 2009 * Modified by William Zhang, 2013 *=================================================================*/ /* Required Packages */ #include #include /* Input Arguments */ #define T_row_IN prhs[0] /* Output Arguments */ #define LD_OUT plhs[0] void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[] ) { /* Declare variables */ double *T_row, *logdet, *r, *lambda, *ghat, *ghatold, *gamma; unsigned int n; int i,j; /* Check for proper number of arguments*/ if (nrhs != 1) { mexErrMsgTxt("1 input arg required."); } else if (nlhs > 2) { mexErrMsgTxt("2 output args required."); } /* Check the dimensions of the arguments */ /* This expects T to be an n+1 by n+1 matrix */ n = mxGetN(T_row_IN)-1; /* Create a matrix for the return argument */ LD_OUT = mxCreateDoubleScalar(0); /* LD_OUT is a pointer to a 1x1 mxArray double initialized to 0. */ /* Assign pointers to some input and output parameters */ logdet = mxGetPr(LD_OUT); T_row = mxGetPr(T_row_IN); /*Allocate space on heap for various variables.*/ r = (double*) mxMalloc(n*sizeof(double)); if (r == NULL) { mexErrMsgTxt("Inadequate space on heap for r."); } lambda = (double*) mxMalloc(n*sizeof(double)); if (lambda == NULL) { mexErrMsgTxt("Inadequate space on heap for lambda."); } ghat = (double*) mxMalloc(n*sizeof(double)); if (ghat == NULL) { mexErrMsgTxt("Inadequate space on heap for ghat."); } ghatold = (double*) mxMalloc(n*sizeof(double)); if (ghatold == NULL) { mexErrMsgTxt("Inadequate space on heap for ghatold."); } gamma = (double*) mxMalloc(n*sizeof(double)); if (gamma == NULL) { mexErrMsgTxt("Inadequate space on heap for gamma."); } /* Define r, the normalized row from T_row(1,2) to T_row(1,n+1) */ for (i = 0; i < n ; i++) { r[i] = T_row[0*(n+1) + (i+1)]/T_row[0*(n+1) + 0]; } /* Initialize the values of lambda and g for the algorithm */ lambda[0] = 1 - pow(r[0],2); ghat[0] = -r[0]; /* Recursion to build g and lambda */ for (i=0 ; i < n-1 ; i++) { /* Calculate gamma */ gamma[i] = -r[i+1]; /* Note that ghat, g_i, etc are i+1 elements long */ for (j=0 ; j < i+1 ; j++) { gamma[i] -= r[j]*ghat[j]; } /* Calculate next ghat */ /* First store ghatold, which is i+1 elts long */ for (j=0 ; j < i+1 ; j++) { ghatold[j] = ghat[j]; } /* The first element of the new ghat is derived from gamma and lambda. */ ghat[0] = gamma[i]/lambda[i]; /* The rest of the new ghat is formed from the old ghat, ghatold, as described on page 599 of the Zohar paper. */ for (j=1 ; j < i+2 ; j++) { ghat[j] = ghatold[j-1] + gamma[i]/lambda[i]*ghatold[i+1-j]; } /* Calculate lambda */ lambda[i+1] = lambda[i] - pow(gamma[i],2)/lambda[i]; } /* There are n lambdas, n g values, and n-1 gammas */ /* Note that n is not the dimension of the matrix, but one less. */ /* This corresponds with Zohar notation. */ /* NOTE ON MEX MATRIX INDEXING */ /* indexing is always (colidx * colLen + rowidx) */ /* Calculate the log determinant for free (essentially) by summing the logs of the lambdas. */ logdet[0] = 0; for (i=0 ; i < n; i++) { logdet[0] += log(lambda[i]); } /* Renormalize based on T_row(1,1) */ logdet[0] += (n+1)*log(T_row[ 0*(n+1) + 0 ]); /* Free allocated arrays */ mxFree(gamma); mxFree(ghatold); mxFree(ghat); mxFree(lambda); mxFree(r); return; } toeblitz-1.0/decomp/logdetToeplitzFastZohar.mexw640000775000442500044270000002200012202534277023435 0ustar cunninghamcunninghamMZ@ !L!This program cannot be run in DOS mode. $vvvv݃v޹vwv܃vvvRichvPEd R"  4p@@']#dP@` .text" `.rdata @@.data`0@.pdata@@@.rsrcP @@.reloc `"@BHHXHhHp WATAUAVAWH)xIHD)@At H  ~ H H fWpHt$PHHH HHD$HHLHHD$XHl$`oLHu H PsHSHD$ Hu H UH5LHu H 9HLHu H HHHD$8Hu H 9E3E3҃|vFIIMI+DD L@AH IA^ED IA^ML AA^ED IA^ML uD;s*MŋK M+A+ADHHA^EAuDsAfA(: THL$ FHDŽ$\ A=fWA$&HHMI+֋M+)$HT$@IAH+LHL$0HT$(H+HDE3BT(3fWAAAIIOI+DDHfDD H IYA\L YI\D YA\L YI\uE;},MAI M+A+BHHYA\uE3E3A|[AAIINI+DDLHDH IHAHDHAHDHAHDHAuE;}!MAK M+A+JHHHAuAyD^+A$;AHD$(IIL$LTAAI+DDH4H I I^+AYB0XD A ^ +AYJ(XL I^+AYB XD A ^ +AYJXL IuD;}GH$MDH+M+I ITE+HHI^+YBBXD Au4+fA()HT$(HL$@^H$HDHHL$0D\t5Ll$XH|$HHt$P($Hl$`HL$ D(D$p($Ht$HHHXuAENfH*HL$8YXCI;I3HL$ )IL$I[0Ik@IsHIA_A^A]A\_% % % % % % % %. @SH  HH HHHuC#H#H H 3H [HHXHhHxL` AUAVAWH 3ML8Q#D>eH%0HXH;t/ 3H4uAtcH  LHH MLHHI;rZH9}tH9EtHMHHEH H HL;uL;tLLIHpHq=SEH=K3eH%0HXH;t3H ut ;>HH uHH u HHH9=t!H tMĺIOH\$@Hl$HH|$PLd$XH A_A^A]HHXHpHxATH0ILXu9u 3ۉXtu7H HtЋ؉D$ LƋI0؉D$ LƋI؉D$ u5u1L3IL3ILMt L3IAӅtu7LƋI#ˋىL$ tH^HtLƋIЋ؉D$ 3ۉ\$ H\$@Ht$HH|$PH0A\H\$Ht$WH IHuoLNjHH\$0Ht$8H _H @SH HH HD$8Hu H~H |HD$8H xjHD$@HTHLD$@HT$8HHL$84HEHL$@"H+WHH [H(GHH(H\$WH HH=HHtHH;rH\$0H _H\$WH HH=HHtHH;rH\$0H _HMZf9t3HcH PAPADDINGXXPADDINGPADDINGXXPADDINGPADDINGXXPADDINGPADDINGXXPADDINGPADDINGXXPAD (toeblitz-1.0/decomp/decompToeplitzFastZohar.mex0000775000442500044270000003245312202534277023102 0ustar cunninghamcunninghamELF4L)4 ($!HH$$QtdRtdttGNUVU`~ڱz6X     ` |CEqX4K _+ "Xt4 , , ?@F    __gmon_start___init_fini__cxa_finalize_Jv_RegisterClassesmexFunctionmexErrMsgTxtmxGetNmxCreateDoubleMatrixmxCreateDoubleScalarmxGetPrmxMalloclogmxFreeliboctinterp.soliboctave.solibcruft.soliblapack.so.3gflibblas.so.3gflibfftw3.so.3libfftw3f.so.3libreadline.so.6libncurses.so.5libdl.so.2libhdf5.so.6libz.so.1libgfortran.so.3libstdc++.so.6libm.so.6libgcc_s.so.1libc.so.6_edata__bss_start_endGLIBC_2.0GLIBC_2.1.3n ii si (           $  US[Lt.=X[ hhhhh h($h0(h8p,h@`0hHPUVSj8u]t4$h<)9s <<9rƃ8[^]US.tt $Ѓ[]Ë$ÐUWVSu à}}x$g$D$EȃEEȉD$$]$F$tEԋF$fE$YtEE$"E2U$ E/M $-}<$E*E$E'Mt11}uܐ4׉19EwUM}uȃ}}ثUEU1ҍEЋ}MljE܋EM1Mf ƃ9~؋M19~M܃1MU1M܋MĉUU2U \9EuȋM܋U؋}̋EE;ULMU}ыM2UvF}Ⱥu}؋}uE)֍u77U9Mu֋u܃}vS}}ȉu܉}؍u܃}؍֋uUU6u9M}7։w͋u܋EUԋMEM}ȋUȋMȋEȃ}}|UȉMفE1UEExMu؍U+Ũ9ЉUu؍H}4֋U؉uu̍u܋u11&M$ݝhE9݅hw֋u܋t$ݝhE1҉UM}Em݅h<$E$4$U$M $ļ[^_]Ã}$y}EEEث $@$h$$n$[UVSZt&Ћu[^]US[ Y[1 input arg required.2 output args required.Inadequate space on heap for r.Inadequate space on heap for lambda.Inadequate space on heap for ghat.Inadequate space on heap for ghatold.Inadequate space on heap for gamma.? ,7DN_nx  o  PL, oooo"2BRbr( GCC: (Ubuntu 4.4.3-4ubuntu5.1) 4.4.3@FmexFunction@  int!95s&ni2 ]1@  1(, A1X 1( K1 W4* Ti4S 4u r4 P4uF44 4u n5Df i6(y j6(==6% $ > $ > : ; $ >  : ; ( .? : ; ' @ : ; I 4: ; I 4: ; I 4: ; I 4: ; I4: ; I I&I /home/cunni/toeplitz/wemmick/toeplitz/decomp/usr/include/octave-3.2.3/octavedecompToeplitzFastZohar.cmex.hmxarray.h@1D ֢~-=k-=k-/k-=k-=m>,5jbNicP(dZH dz?U[9[c^dt<dX)fH9?"}A7 t#9?rU?>|j8T@~~<<~<E~f:0H ;Y-/0:>:0z8$~=.gջ[4444| @FAB LnlhslogdetmxREALnrhsunsigned charshort unsigned intmxArrayplhsghatprhslambdaT_rowmexFunctionlong long unsigned intgammaghatoldmxCOMPLEXlong long intGNU C 4.4.3short intlong double/home/cunni/toeplitz/wemmick/toeplitz/decomp/decompToeplitzFastZohar.cttFuuF''Vu Vu  V Fu ''Fu' 'WuWu W FuBQW-QWV-V!zPzR<R<\P\QWQWPEsPsQ~PP.symtab.strtab.shstrtab.note.gnu.build-id.gnu.hash.dynsym.dynstr.gnu.version.gnu.version_r.rel.dyn.rel.plt.init.text.fini.rodata.eh_frame.ctors.dtors.jcr.dynamic.got.got.plt.data.bss.comment.debug_aranges.debug_pubnames.debug_info.debug_abbrev.debug_line.debug_frame.debug_str.debug_loc$2.o<8  @Ho$Uo@d ,, m LLP v0q|H  H 4( ( , , 0, %Q q " 8"#$4*0$5%  (@.P#4 <3,L     ( ,   (5 K, Z0 h t  7 (   !., :4 ? N, Uj ~"@F  crtstuff.c__CTOR_LIST____DTOR_LIST____JCR_LIST____do_global_dtors_auxcompleted.7021dtor_idx.7023frame_dummy__CTOR_END____FRAME_END____JCR_END____do_global_ctors_auxdecompToeplitzFastZohar.c__i686.get_pc_thunk.bx_GLOBAL_OFFSET_TABLE__DYNAMIC__DTOR_END____dso_handle_finimexErrMsgTxt__bss_start_end__gmon_start___edatamxCreateDoubleMatrix_Jv_RegisterClassesmxMalloc__cxa_finalize@@GLIBC_2.1.3mexFunctionmxGetPr_initlog@@GLIBC_2.0mxGetNmxCreateDoubleScalarmxFreetoeblitz-1.0/decomp/logdetToeplitzFastZohar.mexa640000775000442500044270000001771112202534277023424 0ustar cunninghamcunninghamELF>@@8@      $$Ptd QtdGNUu,L B @  5 hL  S"p}?y3 \__gmon_start____cxa_finalize_Jv_RegisterClassesmexFunctionmexErrMsgTxtmxGetNmxCreateDoubleScalarmxGetPrmxMalloclogmxFreelibmx.solibmex.solibmat.solibm.so.6libstdc++.so.6libpthread.so.0libc.so.6logdetToeplitzFastZohar.mexa64MEXGLIBC_2.2.5 Q ui ui      ( H P X ` h p  x    H:%H52 %4 @%2 h%* h%" h% h% h% h% h% hpHH} HtHÐU= HATSubH=` t H=? jH# L% H L)HHH9s DHH AHz H9rf [A\fH= UHtH HtH= @ÐUHAWAVAUATSHXIH˃tH=~ H=H;EEfWI$HSHEH;GHEEHHEHpIIHu H=yH}PHEHu H=[H}2HHEHu H=:H}IIHu H=H}HEHu H=}}HUHcȃ^A9wAY \HEA fWHEEEDHcLHMLL CfWDBE~8LA Y \HH9uLH<I|HH9uL]A^HEB~.LA^HcAYAXDDHuf(YA^\HEBHD;EHEH}tELuDuLmIA$EHcADjXEA$A9wLuHEEHE9EH*YXEHEH}LH}H}LHX[A\A]A^A_AY \HEA fWHEEAHEHLuDuLmIUHSHH HtH HHHuH[ÐHOH1 input arg required.2 output args required.Inadequate space on heap for r.Inadequate space on heap for lambda.Inadequate space on heap for ghat.Inadequate space on heap for ghatold.Inadequate space on heap for gamma.?;0zRx ,AC P) A    X o0  0 (` oPoooo4o &6FVfvGCC: (GNU) 4.4.6 20120305 (Red Hat 4.4.6-4)GCC: (GNU) 4.4.7 20120313 (Red Hat 4.4.7-3).symtab.strtab.shstrtab.note.gnu.build-id.gnu.hash.dynsym.dynstr.gnu.version.gnu.version_d.gnu.version_r.rela.dyn.rela.plt.init.text.fini.rodata.eh_frame_hdr.eh_frame.ctors.dtors.jcr.data.rel.ro.dynamic.got.got.plt.bss.comment$.o(8 P@00Ho44UoPP8do@s`}(( X X p p   L         0 0X 0X 2 04P (    X p         0    * 8 E [ j x 0     X 0     " ) 2 8@G V jn" \call_gmon_startcrtstuff.c__CTOR_LIST____DTOR_LIST____JCR_LIST____do_global_dtors_auxcompleted.6349dtor_idx.6351frame_dummy__CTOR_END____FRAME_END____JCR_END____do_global_ctors_auxlogdetToeplitzFastZohar.c_fini_GLOBAL_OFFSET_TABLE___dso_handle__DTOR_END____bss_start_end_edata_DYNAMIC_initmxGetPrmxGetN__gmon_start___Jv_RegisterClassesMEXmxCreateDoubleScalar__cxa_finalize@@GLIBC_2.2.5mxMallocmxFreemexFunctionmexErrMsgTxtlog@@GLIBC_2.2.5toeblitz-1.0/decomp/decompToeplitzFastZohar.mexa640000775000442500044270000002146212202534277023413 0ustar cunninghamcunninghamELF>@@8@    HH H $$PtdQtdGNU=qbv(>6 @  5 8L  l"?S 3 __gmon_start____cxa_finalize_Jv_RegisterClassesmexFunctionmexErrMsgTxtmxGetNmxCreateDoubleMatrix_700mxCreateDoubleScalarmxGetPrmxMalloclogmxFreelibmx.solibmex.solibmat.solibm.so.6libstdc++.so.6libpthread.so.0libc.so.6decompToeplitzFastZohar.mexa64MEXGLIBC_2.2.5$Q  ui ui @ @ H P X x             HJH5 % @% h% h% h% h% h% h% h% hp% h`HHM HtHÐU= HATSubH=0 t H= ZH L% H} L)HHH9s DHH] AHR H9r> [A\fH= UHtH HtH= @ÐUHAWAVAUATSHHI΃tH=~ H=I>nAAUUƉHLcfWVI$H;*HI<$HhI>H`EHHEH5IHpHu H=O:H}HEHu H=QH}IIHu H=YH}HEHXHu H=\H}HEHu H=f}}LpH`Hcȃ^A9wHpY 7\LUA  .fWAEEDHcL4HMLLCDfWDHE~8LA AY \HH9uLI<H<HH9uLuA^HUP~-LA^HcAYXADHuf(YA^\HMBHD;MLuH]DmELUIHE^AvEDEHuD)I Ή^A^ӃA9w̓}vN}AuLEI D :A^BAA^Ӄ9wʋEEHHËE艅|@AUUDىMEUDxEDUEEEED؃ADH AEDMf(LUA^ DUA)ACYAA$CY\YBXEÉL L ÉH H Eu;UyEEEDmUU9|vR‰ЋM+MȉM9sӉM<΋MM$DBEDUGD֋}}uAx+MEtIL`+AHA^A9w߃A9vHhH}tKLuDuIELmA$EHcADXEA$A9wLuELhAEL`AnEI*YXEHhH}HX LH}HpHĈ[A\A]A^A_HpY \HM  fWAE}ALXLuH]H]DmLpFUHSHH HtH HHHuH[ÐHH1 input arg required.2 output args required.Inadequate space on heap for r.Inadequate space on heap for lambda.Inadequate space on heap for ghat.Inadequate space on heap for ghatold.Inadequate space on heap for gamma.?; 0zRx ,AC S A @  8 oH  ` `` ooooodoH fvGCC: (GNU) 4.4.6 20120305 (Red Hat 4.4.6-4)GCC: (GNU) 4.4.7 20120313 (Red Hat 4.4.7-3).symtab.strtab.shstrtab.note.gnu.build-id.gnu.hash.dynsym.dynstr.gnu.version.gnu.version_d.gnu.version_r.rela.dyn.rela.plt.init.text.fini.rodata.eh_frame_hdr.eh_frame.ctors.dtors.jcr.data.rel.ro.dynamic.got.got.plt.bss.comment$.o(8 h@HHHoddUo8do@s`}`` 88PPL ( (8 8@ @H HH H` `` 0XX2 @!Hd ` 8 P   ( 8 @ H H `    *( 88 E [ j x   8  P ` @ 0   " )H 2 88@G V jn" call_gmon_startcrtstuff.c__CTOR_LIST____DTOR_LIST____JCR_LIST____do_global_dtors_auxcompleted.6349dtor_idx.6351frame_dummy__CTOR_END____FRAME_END____JCR_END____do_global_ctors_auxdecompToeplitzFastZohar.c_fini_GLOBAL_OFFSET_TABLE___dso_handle__DTOR_END____bss_start_end_edata_DYNAMIC_initmxGetPrmxGetN__gmon_start___Jv_RegisterClassesMEXmxCreateDoubleScalar__cxa_finalize@@GLIBC_2.2.5mxMallocmxFreemexFunctionmexErrMsgTxtmxCreateDoubleMatrix_700log@@GLIBC_2.2.5toeblitz-1.0/decomp/decompToeplitz.m0000775000442500044270000002536412202534277020726 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % John P Cunningham, 2009 % Modified by William Zhang, 2013 % % decompToeplitz() % % Invert a symmetric, real, positive definite Toeplitz matrix % using either inv() or the Trench algorithm, which % uses Zohar 1969, finding the log determinant along the way. % This is slightly different than % Algorithm 4.7.3 of Golub and Van Loan % "Matrix Computations," the 1996 version, because the Zohar % version produces the log determinant essentially for free, which % is often useful in cases where one actually needs the full matrix inverse. % % If Zohar or the MATLAB implementation of the Trench algorithm is called, % the inversion is done in O(n^2) time, which is considerably better % than the O(n^3) that inv() offers. This follows the Trench % algorithm implementation of Zohar 1969. See that paper for % all explanation, as the C code is just an implementation of % the algorithm of p599 of that paper. See Algorithm 4.7.3 of Golub and Van Loan % for an explanation of the MATLAB version of this algorithm, which is slower % due to for loops (slower than inv(), up to about n=300). % Thus, if the C can be used, always use it. % % This function also computes the log determinant, which is % calculated essentially for free as part of the calculation of % the inverse when Zohar is used. This is often useful in applications when one really % needs to represent the inverse of a Toeplitz matrix. If Zohar is not used, the computation % is reasonably costly. % % Usage: [ld, Ti] = decompToeplitz(T, logdetMode, [runMode]); % % Inputs: % T the positive definite symmetric Toeplitz matrix to be inverted % logdetMode if this is set to 1, we skip the building of the inverse for efficiency, % and only compute the log determinant. If it is 0, we still compute both % the log determinant and the inverse % runMode OPTIONAL: to force this inversion to happen using a particular method. This % will NOT roll over to a method (by design), so it will fail if something is amuck. % % Outputs: % ld the log determinant of T, log(det(T)) % Ti the inverted matrix, also symmetric (and persymmetric), NOT TOEPLITZ % % NOTE: We bother with this method because we % need this particular matrix inversion to be % as fast as possible. Thus, no error checking % is done here as that would add needless computation. % Instead, the onus is on the caller to make sure % the matrix is toeplitz and symmetric and real. % % NOTE: Algorithm 4.7.3 in the Golub book has a number of typos % Do not use Alg 4.7.3, use the equations on p197-199. % % Run time tests (testToeplitz) suggest that Zohar is the fastest method across % the board. (2-3 orders of magnitude at sizes from 200-1000). Classical Trench is 13/4n^2 flops, % whereas inv() is probably 2/3n^3. However, the overhead for the MEX calls, etc., % may make the crossover point very MATLAB dependent. % % NOTE: The interested numerical linear algebra geek can dig into all this stuff by % running that testToeplitz to find the crossover % point him/herself. Some tinkering is required. This will run three inversion methods: inv(), decompToeplitzFastZohar (MEX Zohar), % decompToeplitz (Matlab Zohar). This version of decompToeplitz has been pared down to % be more user friendly, so it only includes a call to C/MEX Zohar, inv(), and vectorized MATLAB Zohar. % % MEX NOTE: This function will call a MEX routine. % A try block is included to default back to the best non-MEX solver (either inv() % or a vectorized Trench/Zohar version in native MATLAB), so hopefully the user will % not experience failures. We have also included the compiled .mex in a number of % different architectures (hopefully all). However, the user should really compile % this code on his/her own machine to ensure proper use of MEX. That can be done % from the MATLAB prompt with "mex decompToeplitzFastZohar.c". This C code does nothing % fancy, so if you can run MEX at all, this should work. See also 'help mex' and % 'mexext('all')'. % % MEX NOTE 2: The try-catch block is to hopefully fail MEX gracefully. If the mex % code is not compiled or in the wrong architecture, this will default to the next % best MATLAB implementation, which, depending on n, is inv() or vectorized Trench/Zohar. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ld, Ti] = decompToeplitz(T_row, logdetMode, runMode) n = size(T_row, 2); % Dimension of the matrix % If runMode is specified, just run with that (for testing, etc.). Otherwise, do automated determination of runMode. if nargin < 3 % If runMode is not specified, do the automated determination of runMode, based on runtime tests tryMode = 0; % Always try C/MEX Zohar first. If MEX fails, use inv() or vectorized Zohar Trench in MATLAB. if n < 150 % The n=150 is a rough estimate based on tests on 64 and 32 bit linux systems catchMode = -1; % inv() is the best failure option when n < 150. else catchMode = 1; % n is >= 150, so vectorized Zohar is the best failure option. end else % If specified, force it. tryMode = runMode; catchMode = runMode; end % Invert T with the specified method. try [Ti, ld] = decompToeplitzMode(T_row, logdetMode, tryMode); catch [Ti, ld] = decompToeplitzMode(T_row, logdetMode, catchMode); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Subfunction for try-catch block. This contains the actual Matlab algorithm, and calls the wrapper for the C/MEX function. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Ti, ld] = decompToeplitzMode(T_row, logdetMode, runMode) switch (runMode) % Invert T to produce Ti, the inverse, and ld, the log determinant, by the specified method (or, if logdetMode is on, only the log determinant). case -1 T = convertToeplitz(T_row); % Convert our usual row vector notation into a full Toeplitz matrix. if not(logdetMode) Ti = inv(T); % Just call MATLAB inv(). elseif logdetMode Ti = 0; % Set an arbitrary (incorrect) value for Ti so that everything works smoothly. end ld = logdet(T); % Calculate the log determinant separately. case 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Here is the fast C-MEX implementation of this inversion, using Zohar % invToeplitzFastZohar.c should be MEX-compiled with a command from the MATLAB % prompt of "mex invToeplitzFastZohar.c" This should just work. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if not(logdetMode) [Ti, ld] = decompToeplitzFast(T_row); % Call the MEX wrapper. elseif logdetMode ld = logdetToeplitzFastZohar(T_row); % Call the C/MEX function. Ti = 0; % Set an arbitrary (incorrect) value for Ti so that everything works smoothly. end case 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % A Matlab implementation of the Zohar algorithm. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% scale_factor = T_row(1, 1); % Save the scaling factor. T_row = T_row/scale_factor; % Normalize using the value on the central diagonal so that we can use the Zohar algorithm. n = size(T_row, 2)-1; r = transpose(T_row(1, 2:end)); % Intialize the value of r. In our code we have combined the variables r and rho from the paper, since they contain the same information. lambda = zeros(n, 1); % Initialize the vector of lambda values. % Initialize the algorithm. lambda(1, 1) = 1-r(1, 1)^2; % Obtain the first value of lambda. ghatold(1, 1) = -r(1, 1); % Obtain the first value of ghat, which is the same as g for the length-one case. ghat = ghatold; % Initialize ghat as ghatold. % Recursion to compute the values of gamma and lambda, and g and ghat. In our code we have combined the g and ghat variables, since they contain the same information. for i = 1:n-1 gamma(i, 1) = -(r(i+1, 1) + transpose(r(1:i, 1))*ghatold); % Compute gamma according to the Zohar algorithm. ghat = vertcat(gamma(i, 1)/lambda(i, 1), ghatold + gamma(i, 1)/lambda(i, 1)*ghatold(end:-1:1, 1)); % Compute the next value of ghat according to the Zohar algorithm. lambda(i+1, 1) = lambda(i, 1) - gamma(i, 1)^2/lambda(i, 1); % Compute lambda according to the Zohar algorithm. ghatold = ghat; % Save the new value of ghat in ghatold in preparation for the next recursive step. end if not(logdetMode) % If we're not in logdetMode, then construct and output the inverse, Ti. Ti(1, 1) = 1/lambda(n, 1); % Fill out the first element (upper left) of Ti. Ti(2:n+1, 1) = ghat(n:-1:1, 1)/lambda(n, 1); % Fill out the first column (left) of Ti. Ti(1, 2:n+1) = transpose(Ti(2:n+1, 1)); % Fill out the first row (top) of Ti by symmetry. Ti(n+1, 2:n+1) = Ti(1, n:-1:1); % Fill out the last row (bottom) of Ti by persymmetry and symmetry. Ti(2:n, n+1) = Ti(n:-1:2, 1); % Fill out the last column (right) of Ti by persymmetry and symmetry. % Fill out the rest of Ti using the equations from Zohar and symmetry. for j = 1:floor(n/2); Ti(j+1:n-j+1, j+1) = Ti(j:n-j, j) + 1/lambda(n, 1)*(ghat(n+1-j)*ghat(n+1-j:-1:1+j)-ghat(j:n-j)*ghat(j)); % Apply the Zohar equation to fill in one-fourth of the Ti output matrix. Ti(j+1, j+1:n-j+1) = Ti(j+1:n-j+1, j+1); % Use symmetry of Ti to fill out another quarter of the matrix. Ti(n+1-j:-1:j+1, n+2-(j+1)) = Ti(j+1:n-j+1, j+1); % Use persymmetry to fill out another quarter of the matrix. Ti(n+2-(j+1), n+2-(j+1):-1:n+2-(n-j+1)) = Ti(j+1, j+1:n-j+1); % Use persymmetry to fill out the final quarter of the matrix. end Ti = Ti/scale_factor; % Rescales Ti to produce the correct inverse. elseif logdetMode % If we are in logdetMode, then give an arbitrary filler value for Ti. Ti = 0; % Set an arbitrary (incorrect) value for Ti so that everything works smoothly. end ld = sum(log(lambda)) + log(scale_factor)*(n+1); % Calculates log(det(T)) for free (essentially, when the inverse is also computed) using the lambdas from our recursion. end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%toeblitz-1.0/decomp/decompToeplitzFastZohar.mexmaci640000775000442500044270000002074012202534277024102 0ustar cunninghamcunninghamx__TEXT__text__TEXT.__stubs__TEXT 0 __stub_helper__TEXT ` __const__TEXTP P__cstring__TEXTpp__unwind_info__TEXT\P\__eh_frame__TEXTP__DATA__dyld__DATA__la_symbol_ptr__DATA@H__LINKEDIT  H 8! P  99e3 MA$  0@rpath/libmx.dylib 0@rpath/libmex.dylib 0@rpath/libmat.dylib 84/usr/lib/libstdc++.6.dylib 8/usr/lib/libgcc_s.1.dylib 8{/usr/lib/libSystem.B.dylib&@ )H ASLgAS%_ %^ UHAWAVAUATSHHIt H= | H=H;IL}1DDIfIFI>II~HhH;HEAGHPHHEHHEHu H=5rH}cHEHu H=7TH}EHEHu H=A6H}'HEHu H=KH} HEHu H=UHPt,HEȍ@1@HUD^HUHH9uHEY \HM fWHEHEȃHH@HFHHMHMHIHM1ɺL]DHu fWHuL~uHH}LEf.AY\LMALIHHuۅ~7HH}LEfDfDAIHHuHuLHu^LHu|`HHLELMDfDLUALLUA^LII CY AXA ILIIHuHuLH}TY^\ HHHuHtHHHEȍHHMLHU^A$rOID$HMȍQLEH}@HuB^B^ljAHMHHHRuHPw0HEȍ@HPA$AHEȅM1 HEȍHHUHRt@HLEDfDH}B^ACB^ljAHuHHHuHEȍ@HPA$AăaHHHHGH`HMȍAHEHXHEQAI|1HpHMHufHpH<0D|IALUIHMDfDf(LmL}C^ L}CACYAY\YACXACACC EC C EC HuLUHHIIIsHpHuHXHu|HHHuHHuH9`Hu9fHMȉʉAHu^AHuuH9uHhHHPusfEHEHEȉH*YXMHhH}H}H}xH}oH}HĈ[A\A]A^A_]UHEȉfELufDANXEHhIHOE%L%N%P%R%T%V%X%ZLLL|L pL dLXLLL@?1 input arg required.2 output args required.Inadequate space on heap for r.Inadequate space on heap for lambda.Inadequate space on heap for ghat.Inadequate space on heap for ghatold.Inadequate space on heap for gamma.44 4  XzRx 4 __  ,8D (08@H  .;@Nh~   dyld_stub_binding_helper__dyld_func_lookup_mexFunction_log_mexErrMsgTxt_mxCreateDoubleMatrix_700_mxCreateDoubleScalar_mxFree_mxGetN_mxGetPr_mxMalloctoeblitz-1.0/decomp/logdetToeplitzFastZohar.mexmaci640000775000442500044270000002064012202534277024110 0ustar cunninghamcunninghamx__TEXT__text__TEXT  __stubs__TEXT * __stub_helper__TEXT T __const__TEXTP P__cstring__TEXTpp__unwind_info__TEXT\P\__eh_frame__TEXTP__DATA__dyld__DATA__la_symbol_ptr__DATA8H__LINKEDIT  @ ! P  k'6پI$  0@rpath/libmx.dylib 0@rpath/libmex.dylib 0@rpath/libmat.dylib 84/usr/lib/libstdc++.6.dylib 8/usr/lib/libgcc_s.1.dylib 8{/usr/lib/libSystem.B.dylib&8 )@ ASLAS%%UHAWAVAUATSHHHIt H= | H=%H;jIL}fNIHUHEH;IHEAGHEHH3HEHu H=$HIMu H=HHEHu H=HIMu H=HHHu H=HEt0HE@1HUD^HUHH9uHEY \AfWHEHE[HH@HFHHMHMHIHM1ɺ|HfH} fWL~fHLELMAAY\LIIHu߅~(HMLMАA AIIHuLA^LH}|CHIMLUfLA^LMI CY AX A IIIHuALTY^\A HHH}H|HHHEHHEujfEHEHEH*YXMHEH}LuH}lLdH}HH[A\A]A^A_]MHEAAfEMDAEEXEHEIIZE%>%@%B%D%F%H%JL L L LLLL?1 input arg required.2 output args required.Inadequate space on heap for r.Inadequate space on heap for lambda.Inadequate space on heap for ghat.Inadequate space on heap for ghatold.Inadequate space on heap for gamma. 44 4  XzRx 4`  __ (4@ (08@  $ .0 ;@Ndlt}   dyld_stub_binding_helper__dyld_func_lookup_mexFunction_log_mexErrMsgTxt_mxCreateDoubleScalar_mxFree_mxGetN_mxGetPr_mxMalloctoeblitz-1.0/test/0000775000442500044270000000000012206220476015236 5ustar cunninghamcunninghamtoeblitz-1.0/test/test_solveToeplitz.m0000775000442500044270000000434412202534277021351 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This function takes an n x n Toeplitz matrix T and a nx1 column vector b as input, and returns a structured array with rows containing accuracy results and timings for various methods of calculating solving the system Tx = b for x, and columns corresponding to i, the sizes of the different matrices tested. For each method, the upper row contains the accuracy results for the solution x, while the lower row contains the timings for the solution. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [new_vals_times] = test_solveToeplitz(T_row, b) new_vals_times = zeros(6, 1); % Initialize the vector of results with zeros. % Call solveToeplitz and specify the default Matlab mldivide algorithm. t = cputime; x = solveToeplitz(T_row, b, 'Mldivide'); new_vals_times(2, 1) = cputime-t; % Times the default algorithm. new_vals_times(1, 1) = norm(x(:) - x(:)); % Compares accuracy to built-in Matlab algorithm. % Call solveToeplitz and specify the Levinson-Durbin algorithm. t = cputime; xLD = solveToeplitz(T_row, b, 'LD'); new_vals_times(4, 1) = cputime-t; % Times my algorithm. new_vals_times(3, 1) = norm(x(:) - xLD(:)); % Compares accuracy to built-in Matlab algorithm. % Call solveToeplitz and specify the Conjugate Gradient method. t = cputime; xCG = solveToeplitz(T_row, b, 'PCG_None'); new_vals_times(6, 1) = cputime-t; % Times my algorithm. new_vals_times(5, 1) = norm(x(:) - xCG(:)); % Compares accuracy to built-in Matlab algorithm. % Call solveToeplitz and specify the Preconditioned Conjugate Gradient method with automatic choice of preconditioner. t = cputime; xPCG = solveToeplitz(T_row, b, 'PCG_Auto'); new_vals_times(8, 1) = cputime-t; % Times my algorithm. new_vals_times(7, 1) = norm(x(:) - xPCG(:)); % Compares accuracy to built-in Matlab algorithm. endtoeblitz-1.0/test/test_traceToeplitz.m0000775000442500044270000000416512202534277021320 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This compares the performance of traceToeplitz.m to the built-in trace function for a matrix of size n. For each method, the first row (among the rows its data occupies) contains the norm of the difference between a method and explicitly computing the product and then taking the trace, and the second contains the computational time used. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [new_vals_times] = test_traceToeplitz(T_row) n = size(T_row, 2); new_vals_times = zeros(8, 1); % Initializes new_vals_times. It is faster to modify a matrix than to concatenate two. A = randn(n, n); % This line generates an arbitrary random matrix, A, to use in computing the trace of the product of a Toeplitz matrix with an arbitrary matrix. t = cputime; tr = traceToeplitz(A, T_row, T_row, -1); new_vals_times(2, 1) = cputime-t; % Times the built-in Matlab algorithm. new_vals_times(1, 1) = norm(tr - tr); % Compares accuracy to built-in Matlab algorithm. t = cputime; trT = traceToeplitz(A, T_row, T_row, 0); new_vals_times(4, 1) = cputime-t; % Times my algorithm with automated choice of implementation. new_vals_times(3, 1) = norm(tr - trT); % Compares accuracy to built-in Matlab algorithm. t = cputime; trC = traceToeplitz(A, T_row, T_row, 1); new_vals_times(6, 1) = cputime-t; % Times my algorithm with the fast C implementation. new_vals_times(5, 1) = norm(tr - trC); % Compares accuracy to built-in Matlab algorithm. t = cputime; trM = traceToeplitz(A, T_row, T_row, 2); new_vals_times(8, 1) = cputime-t; % Times my algorithm with the backup Matlab implementation. new_vals_times(7, 1) = norm(tr - trM); % Compares accuracy to built-in Matlab algorithm. endtoeblitz-1.0/test/test_decompToeplitz.m0000775000442500044270000000657012202534277021473 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This function takes a Toeplitz matrix T and its dimension i as input, and returns a structured array with rows containing accuracy results and timings for various methods of calculating the inverse, and columns corresponding to i, the sizes of the different matrices tested. For each method, the upper row contains the accuracy results for the inverse, while the middle row contains the accuracy results for the log determinant, and the lower row contains the timings for the inverse. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[new_vals_times] = test_decompToeplitz(T_row, i) new_vals_times = zeros(12, 1); % Initialize the vector of results with zeros. % Call decompToeplitz(T) and let the code figure out the best method automatically. t0 = clock; [TiInvT, ldT] = decompToeplitz(T_row, 0); new_vals_times(3, 1) = etime(clock,t0); % Saves the time taken to compute the inverse and log determinant in the third row containing information for this method. % Call to decompToeplitz(T,-1) and run the built-in Matlab functions, simple inv() and logdet(). t0 = clock; [TiInv, ld] = decompToeplitz(T_row, 0, -1); tInv(i) = etime(clock,t0); new_vals_times(6, 1) = etime(clock,t0); % Saves the time taken to compute the inverse and log determinant in the third row containing information for this method. % Call to decompToeplitz(T,0) and run the MEX implementation of Zohar. t0 = clock; [TiInvT0, ldT0] = decompToeplitz(T_row, 0, 0); new_vals_times(9, 1) = etime(clock,t0); % Saves the time taken to compute the inverse and log determinant in the third row containing information for this method. % Call to decompToeplitz(T,1) and run the native Matlab implementation of Zohar. t0 = clock; [TiInvT1, ldT1] = decompToeplitz(T_row, 0, 1); new_vals_times(12, 1) = etime(clock,t0); % Saves the time taken to compute the inverse and log determinant in the third row containing information for this method. new_vals_times(1, 1) = norm(TiInv(:) - TiInvT(:)); % Frobenius norm (The (:) linearizes the matrix into a vector and allows computation of the norm.) new_vals_times(4, 1) = norm(TiInv(:) - TiInv(:)); % Frobenius norm (The (:) linearizes the matrix into a vector and allows computation of the norm.) new_vals_times(7, 1) = norm(TiInv(:) - TiInvT0(:)); % Frobenius norm (The (:) linearizes the matrix into a vector and allows computation of the norm.) new_vals_times(10, 1) = norm(TiInv(:) - TiInvT1(:)); % Frobenius norm (The (:) linearizes the matrix into a vector and allows computation of the norm.) new_vals_times(2, 1) = norm(ld - ldT); % Since the ld is a scalar, the norm is much simpler to obtain. new_vals_times(5, 1) = norm(ld - ld); % Since the ld is a scalar, the norm is much simpler to obtain. new_vals_times(8, 1) = norm(ld - ldT0); % Since the ld is a scalar, the norm is much simpler to obtain. new_vals_times(11, 1) = norm(ld - ldT1); % Since the ld is a scalar, the norm is much simpler to obtain. endtoeblitz-1.0/test/test_plot.m0000775000442500044270000001044012202535400017423 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This function plots results, given in data_mat, on a log axis against the vector/matrix sizes, given in sizes_vec, on a linear axis. data_labels contains strings which are used to label to axes. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This main function serves primarily to parse the plot_mode input and use it to call test_plotmode appropriately. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [] = test_plot(full_plot_mode, sizes_vec, data_mat, data_labels, yaxis_label, save_flag, folder, file_name, fig_save) if strcmp(lower(full_plot_mode), 'all') test_plotmode('linear', sizes_vec, data_mat, data_labels, yaxis_label, save_flag, folder, file_name, fig_save); test_plotmode('semilogy', sizes_vec, data_mat, data_labels, yaxis_label, save_flag, folder, file_name, fig_save); test_plotmode('loglog', sizes_vec, data_mat, data_labels, yaxis_label, save_flag, folder, file_name, fig_save); else test_plotmode(lower(full_plot_mode), sizes_vec, data_mat, data_labels, yaxis_label, save_flag, folder, file_name, fig_save); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % test_plotmode actually produces a plot from the input passed to it by test_plot. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [a_graph] = test_plotmode(plot_mode, sizes_vec, data_mat, data_labels, yaxis_label, save_flag, folder, file_name, fig_save) data_series_num = size(data_mat, 1); % Figures out how many data series are being plotted; in this package, it will always be 3 or 4. a_graph = figure; % Creates a figure. % If there are three data series, produce the appropriate figure based on plot_mode. if data_series_num == 3 if strcmp(lower(plot_mode), 'linear') plot(sizes_vec, data_mat(1, 1:end), 'r', sizes_vec, data_mat(2, 1:end), 'b', sizes_vec, data_mat(3, 1:end), 'm', 'LineWidth',3); elseif strcmp(lower(plot_mode), 'semilogy') semilogy(sizes_vec, data_mat(1, 1:end), 'r', sizes_vec, data_mat(2, 1:end), 'b', sizes_vec, data_mat(3, 1:end), 'm', 'LineWidth',3); elseif strcmp(lower(plot_mode), 'loglog') loglog(sizes_vec, data_mat(1, 1:end), 'r', sizes_vec, data_mat(2, 1:end), 'b', sizes_vec, data_mat(3, 1:end), 'm', 'LineWidth',3); end % If there are four data series, produce the appropriate figure based on plot_mode. elseif data_series_num == 4 if strcmp(lower(plot_mode), 'linear') plot(sizes_vec, data_mat(1, 1:end), 'r', sizes_vec, data_mat(2, 1:end), 'b', sizes_vec, data_mat(3, 1:end), 'm', sizes_vec, data_mat(4, 1:end), 'g', 'LineWidth',3); elseif strcmp(lower(plot_mode), 'semilogy') semilogy(sizes_vec, data_mat(1, 1:end), 'r', sizes_vec, data_mat(2, 1:end), 'b', sizes_vec, data_mat(3, 1:end), 'm', sizes_vec, data_mat(4, 1:end), 'g', 'LineWidth',3); elseif strcmp(lower(plot_mode), 'loglog') loglog(sizes_vec, data_mat(1, 1:end), 'r', sizes_vec, data_mat(2, 1:end), 'b', sizes_vec, data_mat(3, 1:end), 'm', sizes_vec, data_mat(4, 1:end), 'g', 'LineWidth',3); end end % Label the axes appropriately and set the font to a fixed size. xlabel('Matrix Size','FontSize',18); ylabel(yaxis_label,'FontSize',18); set(gca,'fontsize',18); % Label the legend appropriately based on the number of data series. if data_series_num == 3 legend(char(data_labels(1, 1)), char(data_labels(2, 1)), char(data_labels(3, 1))); elseif data_series_num == 4 legend(char(data_labels(1, 1)), char(data_labels(2, 1)), char(data_labels(3, 1)), char(data_labels(4, 1))) end % Either saves or displays the plot, based on the value of save_flag. if save_flag print(a_graph, strcat(folder, file_name), '-dpdf'); if fig_save print(a_graph, strcat(folder, file_name), '-dfig'); end close all else a_graph; end end toeblitz-1.0/test/test_solveToeplitzPCG.m0000775000442500044270000000570212202534277021702 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This function takes an n x n Toeplitz matrix T and a nx1 column vector b as input, and returns a structured array with rows containing accuracy results and timings for various preconditioners used for solving the system Tx = b for x by the conjugate gradient method, and columns corresponding to i, the sizes of the different matrices tested. For each method, the upper row contains the accuracy results for the solution x, while the middle row contains the timings for the solution, and the lower row contains the number of iterations required for the solution to converge. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [new_vals_times] = test_solveToeplitzPCG(T_row, b) new_vals_times = zeros(12, 1); % Intialize new_vals_times as a vector of zeros. x = solveToeplitzMldivide(T_row, b); % Get a standard against which to measure the accuracy of our preconditioned conjugate gradient methods. % Call solveToeplitzPCG and specify the unconditioned conjugate gradient method. t = cputime; [xPCG, k] = solveToeplitzPCG(T_row, b, 'None'); new_vals_times(2, 1) = cputime-t; % Times the default algorithm. new_vals_times(1, 1) = norm(x(:) - xPCG(:)); % Compares accuracy to built-in Matlab algorithm. new_vals_times(3, 1) = k; % Stores the number of iterations required to converge. % Call solveToeplitzPCG and specify built-in MATLAB preconditioned conjugate gradient method. t = cputime; [xPCG, k] = solveToeplitzPCG(T_row, b, 'Matlab'); new_vals_times(5, 1) = cputime-t; % Times the default algorithm. new_vals_times(4, 1) = norm(x(:) - xPCG(:)); % Compares accuracy to built-in Matlab algorithm. new_vals_times(6, 1) = k; % Stores the number of iterations required to converge. % Call solveToeplitzPCG and specify the Chan preconditioned conjugate gradient method. t = cputime; [xPCG, k] = solveToeplitzPCG(T_row, b, 'Chan'); new_vals_times(8, 1) = cputime-t; % Times the default algorithm. new_vals_times(7, 1) = norm(x(:) - xPCG(:)); % Compares accuracy to built-in Matlab algorithm. new_vals_times(9, 1) = k; % Stores the number of iterations required to converge. % Call solveToeplitzPCG and specify the Circulant Embedding preconditioned conjugate gradient method. t = cputime; [xPCG, k] = solveToeplitzPCG(T_row, b, 'Circulant_Embedding'); new_vals_times(11, 1) = cputime-t; % Times the default algorithm. new_vals_times(10, 1) = norm(x(:) - xPCG(:)); % Compares accuracy to built-in Matlab algorithm. new_vals_times(12, 1) = k; % Stores the number of iterations required to converge. endtoeblitz-1.0/traceToeplitz.m0000775000442500044270000001205212202534277017274 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This generates the trace of the product of a Toeplitz matrix T, and an % arbitrary square matrix of compatible size, A. A wrapper function tries to use the C implementation, % but falls back on a Matlab implementation if that fails for any reason. % % The calling syntax is % % traceToeplitz(A, T_row, T_col, runMode) % % where % A is the arbitrary matrix, % T_row is the first row of the Toeplitz matrix, % T_col is the (optional) first column of the Toeplitz matrix, written as a row vector, and % runMode is an optional argument, giving the method by which to calculate the trace, with the following accepted values: % % runMode = -1 This calls the default MATLAB matrix multiplication and trace functions. It's very slow and should only be used for testing. % runMode = 0 This is the default method, which first tries to use the C/MEX implementation of our algorithm, but uses the backup MATLAB code if it fails. % runMode = 1 This is for testing only. It forces the code into using the C/MEX implementation of our algorithm with no backup. % runMode = 2 This forces the code into using the MATLAB implementation of our algorithm with no backup. % % traceToeplitz has a single output argument, % % Tr, which gives the desired trace. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This function mainly serves to parse the input. A, the arbitrary matrix, and T_row, the first row of our Toeplitz matrix, are both required. However, T_col, the first column of our Toeplitz matrix, is not required. When it is absent, we simply assume that the matrix is symmetric. In addition, runMode, which specifies which method to use to compute the trace, is optional. It defaults to case 0, which attempts to use fast C code, and resorts to a Matlab backup if the C code fails. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ Tr ] = traceToeplitz(A, T_row, T_col_or_perhaps_runMode, runMode) if nargin < 2 % In this case, we have too few input arguments, and at least one required argument is missing. error('Not enough input arguments') elseif nargin == 2 % In this case, both T_col and runMode are empty. So we assume the Toeplitz matrix is symmetric and use the default case 0 runMode. T_col = T_row; runMode = 0; elseif nargin == 3 % In this case, we need to determine whether the optional argument given is a runMode or a T_col. if size(T_row) == [1 1] % If we are looking at trivial 1x1 matrices, then it is difficult to distinguish T_col from runMode, and therefore we do not allow omitted arguments in this case. error('Not enough input arguments. For trivial 1x1 matrices, T_row, T_col, and runMode must all be explicitly specified to prevent confusion.') elseif size(T_row) == size(T_col_or_perhaps_runMode) % If T_row is the same size as the given optional argument, we can safely assume that it is the first column of our Toeplitz matrix. T_col = T_col_or_perhaps_runMode; runMode = 0; else % Otherwise, the third argument given must be the runMode. T_col = T_row; runMode = T_col_or_perhaps_runMode; end elseif nargin == 4; T_col = T_col_or_perhaps_runMode; end Tr = traceToeplitzMode(A, T_row, T_col, runMode); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This function calls the actual code which will do computations based on the runMode passed to it by traceToeplitz. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ Tr ] = traceToeplitzMode(A, T_row, T_col, runMode) switch(runMode) case -1 % Trivial case: this uses built-in Matlab functions to compute the trace of the product. This is intended primarily for testing. T = toeplitz(T_col, T_row); Tr = trace(A*T); case 0 % Attempts to use the faster C code, but defaults back to the backup Matlab code if anything goes wrong. This is the default case. try % Here it attempts to call the MEX function compiled from the C source code. Tr = traceToeplitzFast(A, T_row, T_col); catch % Here it falls back on the native Matlab code when the C code fails. Tr = traceToeplitzBackup(A, T_row, T_col); end case 1 % Forces the use of the C code. If it fails, then it really fails. This mode is intended mostly for testing. Tr = traceToeplitzFast(A, T_row, T_col); case 2 % Forces the use of the Matlab code. This is more reliable across systems since Matlab is an interpreted language, but is slower. This is intended mostly for testing. Tr = traceToeplitzBackup(A, T_row, T_col); end endtoeblitz-1.0/solve/0000775000442500044270000000000012206220476015407 5ustar cunninghamcunninghamtoeblitz-1.0/solve/solveToeplitzPCG.m0000775000442500044270000000616312202534277021016 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This solves a system Tx = b for x, where x and b are column vectors and T is a positive definite toeplitz matrix. It uses the Conjugate Gradient method from page 529 of "Matrix Computations" by Golub. It also has a second (optional) output argument, k, which tells us the number of iterations it took to converge to our result. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x, k] = solveToeplitzPCG(T_row, b, preconditioner) n = size(T_row, 2); % Get the size of the matrix we are dealing with. epsilon = 10^(-6); % Set a reasonable error bound for use in stopping condition. kmax = 400; % Set a maximum number of iterations. If it's taking more than 400, something is probably wrong. % Here we automatically set the preconditioner if not specified. if nargin < 3 preconditioner = 'auto'; end % Here we define our preconditioner, m. m is actually just the first column of our full circulant preconditioner M, but for speed we omit the rest of the matrix from explicit representation. if strcmp(lower(preconditioner), 'circulant_embedding') m = PCG_Circulant_Embedding(T_row); elseif or(strcmp(lower(preconditioner), 'chan'), strcmp(lower(preconditioner), 'auto')) m = PCG_Chan(T_row); elseif strcmp(lower(preconditioner), 'none') m = zeros(n, 1); m(1, 1) = 1; end if strcmp(lower(preconditioner), 'matlab') [x, foo1, foo2, k] = pcg(toeplitz(T_row), b); % Calls Matlab's default preconditioned conjugate gradient algorithm. We include this here primarily for testing purposes. Note foo used instead of ~ for octave compatibility else x = zeros(n, 1); % Here we choose an initial guess for x. The solution in this method is fairly insensitive to our initial guess, so we just choose the zero vector. k = 0; % Here we initialize k, the number of iterations, at zero. % Initialize the values of a bunch of other variables. r = b - multiplyToeplitz(T_row, x); d = multiplyinvCirculant(m, r); delta_new = transpose(r)*d; delta_old = delta_new; %The following loop runs the preconditioned conjugate gradient method. while and(sqrt(delta_old) > epsilon, k < kmax) q = multiplyToeplitz(T_row, d); alpha = delta_new/(transpose(d)*q); x = x + alpha*d; if mod(k, 50) == 0 r = b - multiplyToeplitz(T_row, x); else r = r - alpha*q; end s = multiplyinvCirculant(m, r); delta_old = delta_new; delta_new = transpose(r)*s; beta = delta_new/delta_old; d = s + beta*d; k = k + 1; end % Prints a warning message if the maximum number of iterations is reached. if k == kmax fprintf('WARNING: %d (maximum number) iterations taken, the residual is %g.', k, delta_old); end end end toeblitz-1.0/solve/PCG_Circulant_Embedding.m0000775000442500044270000000221712202534277022150 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % Given a symmetric positive definite Toeplitz matrix, T, this function returns the first column of the circulant preconditioner formed by circulant embedding as in 2.1.3 (Preconditioners by embedding) of Conjugate gradient methods for Toeplitz systems (1996) by Raymond H. Chan and Michael K. Ng. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [r] = PCG_Circulant_Embedding(T_row) b = vertcat(T_row(1, 1), transpose(T_row(1, end:-1:2))); % Generate the first column of the matrix B, which is embedding in a circulant matrix along with our original matrix T. r = b + transpose(T_row(1, 1:end)); % Add the first column of T to the first column of B to generate the first column of our desired preconditioner. endtoeblitz-1.0/solve/solveToeplitzLD.m0000775000442500044270000000347612202534277020710 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This solves a system Tx = b for x, where x and b are column vectors and T is a positive definite Toeplitz matrix. It uses the Levinson-Durbin algorithm. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x] = solveToeplitzLD(T_row, b) f_v = 1/T_row(1, 1); % This is the first forward vector. b_v = 1/T_row(1, 1); % This is the first backward vector. x = b(1, 1)/T_row(1, 1); % This is the first vector in the sequence x_n leading up to x_N = x, x_1. for i = 2:size(T_row, 2) %Obtaining the forward and backward vectors. f_error = T_row(1, i:-1:2)*f_v; %Computing the error terms that arise when we multiply T_n by the temporary vector. b_error = T_row(1, 2:i)*b_v; x_error = T_row(1, i:-1:2)*x; f_v = (1/(1-f_error*b_error))*(vertcat(f_v, 0) - f_error*vertcat(0, b_v)); %Computes new values for the forward and backwards vectors, as well as the next vector in the x_n sequence. b_v = f_v(end:-1:1); %b_v = (1/(1-f_error*b_error))*b_temp - (b_error/(1-f_error*b_error))*f_temp; %This is commented out, since we are working with symmetric matrices; so we can obtain the backward vector simply by reversing the forward vector. x = vertcat(x, 0) + (b(i, 1)-x_error)*b_v; %Ti*f_v % You can use this to check if the forward vectors are correct; this line of code should produce a vector with a one in the first entry and all zeroes following if they are. end endtoeblitz-1.0/solve/PCG_Chan.m0000775000442500044270000000323112202534277017134 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This function builds the first column of the Chan circulant preconditioner as described in 2.1.2 (T. Chan's Preconditioner) of Conjugate gradient methods for Toeplitz systems (1996) by Raymond H. Chan and Michael K. Ng. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[c] = PCG_Chan(T_row) n = size(T_row, 2); % Gets the size of the Toeplitz matrix we are dealing with. c(1, 1) = T_row(1, 1); % Obtains the first entry of our column vector directly from the Toeplitz matrix. % Non-vectorized version. %for i = 1:n-1 % Loops over all following entries. % c(i+1, 1) = ((n-i)*T_row(1, i+1) + i*T_row(1, n-i+1))/n; % Computes following entries according to equation 2.4 from Chan (1996). %end % Vectorized version. This should be somewhat faster, but is more costly in terms of memory requirements. n_vec = transpose(zeros(1, n-1)); % Initializes n_vec as a vector of all zeros. n_vec(1:n-1, 1) = n; % Sets the value of every entry of n_vec to n. i_vec = transpose([1:n-1]); % Creates a vector with the desired indices we want to loop over. c(2:n, 1) = ((n_vec-i_vec).*transpose(T_row(1, 2:n)) + i_vec.*transpose(T_row(1, n:-1:2)))/n; % % Computes following entries according to equation 2.4 from Chan (1996). endtoeblitz-1.0/solve/solveToeplitzMldivide.m0000775000442500044270000000130612202534277022134 0ustar cunninghamcunningham%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Toeblitz: A Fast Toolkit for Manipulating Toeplitz Matrices % Copyright (C) 2013 William B. Zhang and John P. Cunningham (see full notice in README) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % William Zhang, 2013 % This solves a system Tx = b for x, where x and b are column vectors and T is a toeplitz matrix. It uses the built-in Matlab method. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x] = solveToeplitzMldivide(T_row, b) x = mldivide(toeplitz(T_row), b); endtoeblitz-1.0/license.txt0000775000442500044270000007660212202534277016463 0ustar cunninghamcunninghamGNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007 Copyright © 2007 Free Software Foundation, Inc. Everyone is permitted to copy and distribute verbatim copies of this license document, but changing it is not allowed. Preamble The GNU General Public License is a free, copyleft license for software and other kinds of works. The licenses for most software and other practical works are designed to take away your freedom to share and change the works. By contrast, the GNU General Public License is intended to guarantee your freedom to share and change all versions of a program--to make sure it remains free software for all its users. We, the Free Software Foundation, use the GNU General Public License for most of our software; it applies also to any other work released this way by its authors. You can apply it to your programs, too. When we speak of free software, we are referring to freedom, not price. Our General Public Licenses are designed to make sure that you have the freedom to distribute copies of free software (and charge for them if you wish), that you receive source code or can get it if you want it, that you can change the software or use pieces of it in new free programs, and that you know you can do these things. To protect your rights, we need to prevent others from denying you these rights or asking you to surrender the rights. Therefore, you have certain responsibilities if you distribute copies of the software, or if you modify it: responsibilities to respect the freedom of others. For example, if you distribute copies of such a program, whether gratis or for a fee, you must pass on to the recipients the same freedoms that you received. You must make sure that they, too, receive or can get the source code. And you must show them these terms so they know their rights. Developers that use the GNU GPL protect your rights with two steps: (1) assert copyright on the software, and (2) offer you this License giving you legal permission to copy, distribute and/or modify it. For the developers' and authors' protection, the GPL clearly explains that there is no warranty for this free software. For both users' and authors' sake, the GPL requires that modified versions be marked as changed, so that their problems will not be attributed erroneously to authors of previous versions. Some devices are designed to deny users access to install or run modified versions of the software inside them, although the manufacturer can do so. This is fundamentally incompatible with the aim of protecting users' freedom to change the software. The systematic pattern of such abuse occurs in the area of products for individuals to use, which is precisely where it is most unacceptable. Therefore, we have designed this version of the GPL to prohibit the practice for those products. If such problems arise substantially in other domains, we stand ready to extend this provision to those domains in future versions of the GPL, as needed to protect the freedom of users. Finally, every program is threatened constantly by software patents. States should not allow patents to restrict development and use of software on general-purpose computers, but in those that do, we wish to avoid the special danger that patents applied to a free program could make it effectively proprietary. To prevent this, the GPL assures that patents cannot be used to render the program non-free. The precise terms and conditions for copying, distribution and modification follow. TERMS AND CONDITIONS 0. Definitions. “This License” refers to version 3 of the GNU General Public License. “Copyright” also means copyright-like laws that apply to other kinds of works, such as semiconductor masks. “The Program” refers to any copyrightable work licensed under this License. Each licensee is addressed as “you”. “Licensees” and “recipients” may be individuals or organizations. To “modify” a work means to copy from or adapt all or part of the work in a fashion requiring copyright permission, other than the making of an exact copy. The resulting work is called a “modified version” of the earlier work or a work “based on” the earlier work. A “covered work” means either the unmodified Program or a work based on the Program. To “propagate” a work means to do anything with it that, without permission, would make you directly or secondarily liable for infringement under applicable copyright law, except executing it on a computer or modifying a private copy. Propagation includes copying, distribution (with or without modification), making available to the public, and in some countries other activities as well. To “convey” a work means any kind of propagation that enables other parties to make or receive copies. Mere interaction with a user through a computer network, with no transfer of a copy, is not conveying. An interactive user interface displays “Appropriate Legal Notices” to the extent that it includes a convenient and prominently visible feature that (1) displays an appropriate copyright notice, and (2) tells the user that there is no warranty for the work (except to the extent that warranties are provided), that licensees may convey the work under this License, and how to view a copy of this License. If the interface presents a list of user commands or options, such as a menu, a prominent item in the list meets this criterion. 1. Source Code. The “source code” for a work means the preferred form of the work for making modifications to it. “Object code” means any non-source form of a work. A “Standard Interface” means an interface that either is an official standard defined by a recognized standards body, or, in the case of interfaces specified for a particular programming language, one that is widely used among developers working in that language. The “System Libraries” of an executable work include anything, other than the work as a whole, that (a) is included in the normal form of packaging a Major Component, but which is not part of that Major Component, and (b) serves only to enable use of the work with that Major Component, or to implement a Standard Interface for which an implementation is available to the public in source code form. A “Major Component”, in this context, means a major essential component (kernel, window system, and so on) of the specific operating system (if any) on which the executable work runs, or a compiler used to produce the work, or an object code interpreter used to run it. The “Corresponding Source” for a work in object code form means all the source code needed to generate, install, and (for an executable work) run the object code and to modify the work, including scripts to control those activities. However, it does not include the work's System Libraries, or general-purpose tools or generally available free programs which are used unmodified in performing those activities but which are not part of the work. For example, Corresponding Source includes interface definition files associated with source files for the work, and the source code for shared libraries and dynamically linked subprograms that the work is specifically designed to require, such as by intimate data communication or control flow between those subprograms and other parts of the work. The Corresponding Source need not include anything that users can regenerate automatically from other parts of the Corresponding Source. The Corresponding Source for a work in source code form is that same work. 2. Basic Permissions. All rights granted under this License are granted for the term of copyright on the Program, and are irrevocable provided the stated conditions are met. This License explicitly affirms your unlimited permission to run the unmodified Program. The output from running a covered work is covered by this License only if the output, given its content, constitutes a covered work. This License acknowledges your rights of fair use or other equivalent, as provided by copyright law. You may make, run and propagate covered works that you do not convey, without conditions so long as your license otherwise remains in force. You may convey covered works to others for the sole purpose of having them make modifications exclusively for you, or provide you with facilities for running those works, provided that you comply with the terms of this License in conveying all material for which you do not control copyright. Those thus making or running the covered works for you must do so exclusively on your behalf, under your direction and control, on terms that prohibit them from making any copies of your copyrighted material outside their relationship with you. Conveying under any other circumstances is permitted solely under the conditions stated below. Sublicensing is not allowed; section 10 makes it unnecessary. 3. Protecting Users' Legal Rights From Anti-Circumvention Law. No covered work shall be deemed part of an effective technological measure under any applicable law fulfilling obligations under article 11 of the WIPO copyright treaty adopted on 20 December 1996, or similar laws prohibiting or restricting circumvention of such measures. When you convey a covered work, you waive any legal power to forbid circumvention of technological measures to the extent such circumvention is effected by exercising rights under this License with respect to the covered work, and you disclaim any intention to limit operation or modification of the work as a means of enforcing, against the work's users, your or third parties' legal rights to forbid circumvention of technological measures. 4. Conveying Verbatim Copies. You may convey verbatim copies of the Program's source code as you receive it, in any medium, provided that you conspicuously and appropriately publish on each copy an appropriate copyright notice; keep intact all notices stating that this License and any non-permissive terms added in accord with section 7 apply to the code; keep intact all notices of the absence of any warranty; and give all recipients a copy of this License along with the Program. You may charge any price or no price for each copy that you convey, and you may offer support or warranty protection for a fee. 5. Conveying Modified Source Versions. You may convey a work based on the Program, or the modifications to produce it from the Program, in the form of source code under the terms of section 4, provided that you also meet all of these conditions: a) The work must carry prominent notices stating that you modified it, and giving a relevant date. b) The work must carry prominent notices stating that it is released under this License and any conditions added under section 7. This requirement modifies the requirement in section 4 to “keep intact all notices”. c) You must license the entire work, as a whole, under this License to anyone who comes into possession of a copy. This License will therefore apply, along with any applicable section 7 additional terms, to the whole of the work, and all its parts, regardless of how they are packaged. This License gives no permission to license the work in any other way, but it does not invalidate such permission if you have separately received it. d) If the work has interactive user interfaces, each must display Appropriate Legal Notices; however, if the Program has interactive interfaces that do not display Appropriate Legal Notices, your work need not make them do so. A compilation of a covered work with other separate and independent works, which are not by their nature extensions of the covered work, and which are not combined with it such as to form a larger program, in or on a volume of a storage or distribution medium, is called an “aggregate” if the compilation and its resulting copyright are not used to limit the access or legal rights of the compilation's users beyond what the individual works permit. Inclusion of a covered work in an aggregate does not cause this License to apply to the other parts of the aggregate. 6. Conveying Non-Source Forms. You may convey a covered work in object code form under the terms of sections 4 and 5, provided that you also convey the machine-readable Corresponding Source under the terms of this License, in one of these ways: a) Convey the object code in, or embodied in, a physical product (including a physical distribution medium), accompanied by the Corresponding Source fixed on a durable physical medium customarily used for software interchange. b) Convey the object code in, or embodied in, a physical product (including a physical distribution medium), accompanied by a written offer, valid for at least three years and valid for as long as you offer spare parts or customer support for that product model, to give anyone who possesses the object code either (1) a copy of the Corresponding Source for all the software in the product that is covered by this License, on a durable physical medium customarily used for software interchange, for a price no more than your reasonable cost of physically performing this conveying of source, or (2) access to copy the Corresponding Source from a network server at no charge. c) Convey individual copies of the object code with a copy of the written offer to provide the Corresponding Source. This alternative is allowed only occasionally and noncommercially, and only if you received the object code with such an offer, in accord with subsection 6b. d) Convey the object code by offering access from a designated place (gratis or for a charge), and offer equivalent access to the Corresponding Source in the same way through the same place at no further charge. You need not require recipients to copy the Corresponding Source along with the object code. If the place to copy the object code is a network server, the Corresponding Source may be on a different server (operated by you or a third party) that supports equivalent copying facilities, provided you maintain clear directions next to the object code saying where to find the Corresponding Source. Regardless of what server hosts the Corresponding Source, you remain obligated to ensure that it is available for as long as needed to satisfy these requirements. e) Convey the object code using peer-to-peer transmission, provided you inform other peers where the object code and Corresponding Source of the work are being offered to the general public at no charge under subsection 6d. A separable portion of the object code, whose source code is excluded from the Corresponding Source as a System Library, need not be included in conveying the object code work. A “User Product” is either (1) a “consumer product”, which means any tangible personal property which is normally used for personal, family, or household purposes, or (2) anything designed or sold for incorporation into a dwelling. In determining whether a product is a consumer product, doubtful cases shall be resolved in favor of coverage. For a particular product received by a particular user, “normally used” refers to a typical or common use of that class of product, regardless of the status of the particular user or of the way in which the particular user actually uses, or expects or is expected to use, the product. A product is a consumer product regardless of whether the product has substantial commercial, industrial or non-consumer uses, unless such uses represent the only significant mode of use of the product. “Installation Information” for a User Product means any methods, procedures, authorization keys, or other information required to install and execute modified versions of a covered work in that User Product from a modified version of its Corresponding Source. The information must suffice to ensure that the continued functioning of the modified object code is in no case prevented or interfered with solely because modification has been made. If you convey an object code work under this section in, or with, or specifically for use in, a User Product, and the conveying occurs as part of a transaction in which the right of possession and use of the User Product is transferred to the recipient in perpetuity or for a fixed term (regardless of how the transaction is characterized), the Corresponding Source conveyed under this section must be accompanied by the Installation Information. But this requirement does not apply if neither you nor any third party retains the ability to install modified object code on the User Product (for example, the work has been installed in ROM). The requirement to provide Installation Information does not include a requirement to continue to provide support service, warranty, or updates for a work that has been modified or installed by the recipient, or for the User Product in which it has been modified or installed. Access to a network may be denied when the modification itself materially and adversely affects the operation of the network or violates the rules and protocols for communication across the network. Corresponding Source conveyed, and Installation Information provided, in accord with this section must be in a format that is publicly documented (and with an implementation available to the public in source code form), and must require no special password or key for unpacking, reading or copying. 7. Additional Terms. “Additional permissions” are terms that supplement the terms of this License by making exceptions from one or more of its conditions. Additional permissions that are applicable to the entire Program shall be treated as though they were included in this License, to the extent that they are valid under applicable law. If additional permissions apply only to part of the Program, that part may be used separately under those permissions, but the entire Program remains governed by this License without regard to the additional permissions. When you convey a copy of a covered work, you may at your option remove any additional permissions from that copy, or from any part of it. (Additional permissions may be written to require their own removal in certain cases when you modify the work.) You may place additional permissions on material, added by you to a covered work, for which you have or can give appropriate copyright permission. Notwithstanding any other provision of this License, for material you add to a covered work, you may (if authorized by the copyright holders of that material) supplement the terms of this License with terms: a) Disclaiming warranty or limiting liability differently from the terms of sections 15 and 16 of this License; or b) Requiring preservation of specified reasonable legal notices or author attributions in that material or in the Appropriate Legal Notices displayed by works containing it; or c) Prohibiting misrepresentation of the origin of that material, or requiring that modified versions of such material be marked in reasonable ways as different from the original version; or d) Limiting the use for publicity purposes of names of licensors or authors of the material; or e) Declining to grant rights under trademark law for use of some trade names, trademarks, or service marks; or f) Requiring indemnification of licensors and authors of that material by anyone who conveys the material (or modified versions of it) with contractual assumptions of liability to the recipient, for any liability that these contractual assumptions directly impose on those licensors and authors. All other non-permissive additional terms are considered “further restrictions” within the meaning of section 10. If the Program as you received it, or any part of it, contains a notice stating that it is governed by this License along with a term that is a further restriction, you may remove that term. If a license document contains a further restriction but permits relicensing or conveying under this License, you may add to a covered work material governed by the terms of that license document, provided that the further restriction does not survive such relicensing or conveying. If you add terms to a covered work in accord with this section, you must place, in the relevant source files, a statement of the additional terms that apply to those files, or a notice indicating where to find the applicable terms. Additional terms, permissive or non-permissive, may be stated in the form of a separately written license, or stated as exceptions; the above requirements apply either way. 8. Termination. You may not propagate or modify a covered work except as expressly provided under this License. Any attempt otherwise to propagate or modify it is void, and will automatically terminate your rights under this License (including any patent licenses granted under the third paragraph of section 11). However, if you cease all violation of this License, then your license from a particular copyright holder is reinstated (a) provisionally, unless and until the copyright holder explicitly and finally terminates your license, and (b) permanently, if the copyright holder fails to notify you of the violation by some reasonable means prior to 60 days after the cessation. Moreover, your license from a particular copyright holder is reinstated permanently if the copyright holder notifies you of the violation by some reasonable means, this is the first time you have received notice of violation of this License (for any work) from that copyright holder, and you cure the violation prior to 30 days after your receipt of the notice. Termination of your rights under this section does not terminate the licenses of parties who have received copies or rights from you under this License. If your rights have been terminated and not permanently reinstated, you do not qualify to receive new licenses for the same material under section 10. 9. Acceptance Not Required for Having Copies. You are not required to accept this License in order to receive or run a copy of the Program. Ancillary propagation of a covered work occurring solely as a consequence of using peer-to-peer transmission to receive a copy likewise does not require acceptance. However, nothing other than this License grants you permission to propagate or modify any covered work. These actions infringe copyright if you do not accept this License. Therefore, by modifying or propagating a covered work, you indicate your acceptance of this License to do so. 10. Automatic Licensing of Downstream Recipients. Each time you convey a covered work, the recipient automatically receives a license from the original licensors, to run, modify and propagate that work, subject to this License. You are not responsible for enforcing compliance by third parties with this License. An “entity transaction” is a transaction transferring control of an organization, or substantially all assets of one, or subdividing an organization, or merging organizations. If propagation of a covered work results from an entity transaction, each party to that transaction who receives a copy of the work also receives whatever licenses to the work the party's predecessor in interest had or could give under the previous paragraph, plus a right to possession of the Corresponding Source of the work from the predecessor in interest, if the predecessor has it or can get it with reasonable efforts. You may not impose any further restrictions on the exercise of the rights granted or affirmed under this License. For example, you may not impose a license fee, royalty, or other charge for exercise of rights granted under this License, and you may not initiate litigation (including a cross-claim or counterclaim in a lawsuit) alleging that any patent claim is infringed by making, using, selling, offering for sale, or importing the Program or any portion of it. 11. Patents. A “contributor” is a copyright holder who authorizes use under this License of the Program or a work on which the Program is based. The work thus licensed is called the contributor's “contributor version”. A contributor's “essential patent claims” are all patent claims owned or controlled by the contributor, whether already acquired or hereafter acquired, that would be infringed by some manner, permitted by this License, of making, using, or selling its contributor version, but do not include claims that would be infringed only as a consequence of further modification of the contributor version. For purposes of this definition, “control” includes the right to grant patent sublicenses in a manner consistent with the requirements of this License. Each contributor grants you a non-exclusive, worldwide, royalty-free patent license under the contributor's essential patent claims, to make, use, sell, offer for sale, import and otherwise run, modify and propagate the contents of its contributor version. In the following three paragraphs, a “patent license” is any express agreement or commitment, however denominated, not to enforce a patent (such as an express permission to practice a patent or covenant not to sue for patent infringement). To “grant” such a patent license to a party means to make such an agreement or commitment not to enforce a patent against the party. If you convey a covered work, knowingly relying on a patent license, and the Corresponding Source of the work is not available for anyone to copy, free of charge and under the terms of this License, through a publicly available network server or other readily accessible means, then you must either (1) cause the Corresponding Source to be so available, or (2) arrange to deprive yourself of the benefit of the patent license for this particular work, or (3) arrange, in a manner consistent with the requirements of this License, to extend the patent license to downstream recipients. “Knowingly relying” means you have actual knowledge that, but for the patent license, your conveying the covered work in a country, or your recipient's use of the covered work in a country, would infringe one or more identifiable patents in that country that you have reason to believe are valid. If, pursuant to or in connection with a single transaction or arrangement, you convey, or propagate by procuring conveyance of, a covered work, and grant a patent license to some of the parties receiving the covered work authorizing them to use, propagate, modify or convey a specific copy of the covered work, then the patent license you grant is automatically extended to all recipients of the covered work and works based on it. A patent license is “discriminatory” if it does not include within the scope of its coverage, prohibits the exercise of, or is conditioned on the non-exercise of one or more of the rights that are specifically granted under this License. You may not convey a covered work if you are a party to an arrangement with a third party that is in the business of distributing software, under which you make payment to the third party based on the extent of your activity of conveying the work, and under which the third party grants, to any of the parties who would receive the covered work from you, a discriminatory patent license (a) in connection with copies of the covered work conveyed by you (or copies made from those copies), or (b) primarily for and in connection with specific products or compilations that contain the covered work, unless you entered into that arrangement, or that patent license was granted, prior to 28 March 2007. Nothing in this License shall be construed as excluding or limiting any implied license or other defenses to infringement that may otherwise be available to you under applicable patent law. 12. No Surrender of Others' Freedom. If conditions are imposed on you (whether by court order, agreement or otherwise) that contradict the conditions of this License, they do not excuse you from the conditions of this License. If you cannot convey a covered work so as to satisfy simultaneously your obligations under this License and any other pertinent obligations, then as a consequence you may not convey it at all. For example, if you agree to terms that obligate you to collect a royalty for further conveying from those to whom you convey the Program, the only way you could satisfy both those terms and this License would be to refrain entirely from conveying the Program. 13. Use with the GNU Affero General Public License. Notwithstanding any other provision of this License, you have permission to link or combine any covered work with a work licensed under version 3 of the GNU Affero General Public License into a single combined work, and to convey the resulting work. The terms of this License will continue to apply to the part which is the covered work, but the special requirements of the GNU Affero General Public License, section 13, concerning interaction through a network will apply to the combination as such. 14. Revised Versions of this License. The Free Software Foundation may publish revised and/or new versions of the GNU General Public License from time to time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or concerns. Each version is given a distinguishing version number. If the Program specifies that a certain numbered version of the GNU General Public License “or any later version” applies to it, you have the option of following the terms and conditions either of that numbered version or of any later version published by the Free Software Foundation. If the Program does not specify a version number of the GNU General Public License, you may choose any version ever published by the Free Software Foundation. If the Program specifies that a proxy can decide which future versions of the GNU General Public License can be used, that proxy's public statement of acceptance of a version permanently authorizes you to choose that version for the Program. Later license versions may give you additional or different permissions. However, no additional obligations are imposed on any author or copyright holder as a result of your choosing to follow a later version. 15. Disclaimer of Warranty. THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM “AS IS” WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. 16. Limitation of Liability. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES. 17. Interpretation of Sections 15 and 16. If the disclaimer of warranty and limitation of liability provided above cannot be given local legal effect according to their terms, reviewing courts shall apply local law that most closely approximates an absolute waiver of all civil liability in connection with the Program, unless a warranty or assumption of liability accompanies a copy of the Program in return for a fee. END OF TERMS AND CONDITIONS