Sparse Matrices Basic Support - II
1. Sparse-specific function: sparse(), nnz(), find().
2. Indexing & basic manipulation functions - subsref, subsasgn, size, numel, find, etc.
3. Element-wise operations (arithmetic and functions).
4. Matrix operations (matrix-scalar, matrix-vector, matrix-matrix).
spdiags with low number of digits
In a larger code I want to construct an mp sparse matrix in two steps - first the off-diagonal entries and then add the diagonal.
However, using mp.Digits(1) the code breaks down with the error message
" Error using mp/sparse
Index exceeds matrix dimensions.
Error in mpspdiags
Error in mp/spdiags (line 4541)
[varargout{1:nargout}] = mpspdiags(varargin{:}); "
I tried to simplify the matrix creation process and got to the following:
mp.Digits(1);
x = zeros(n,'mp');
identity = spdiags(ones(n,1,'mp'),0,x);
This runs fine for n=10, 100 but breaks down for n=101 and many others (with the error message above). Taking mp.Digits(2); and n=1311 or n = 2520 reproduces the same error message. So does mp.Digits(3); and n=4140.
super quick reversion to double
Is there a quick & simple but revertible way to convert an extensive Matlab program, developed in multiprecision, back to double precision? For example, mp.Digits(16) ?
This question arises because code is sometimes developed in multiprecision for the purpose of eliminating quantization artifacts (finite precision) as a possible error source. A natural question arises as to the exact tradeoff of speed versus precision. This information is required for coherently answering one's boss during an entire development phase, preparing talks & papers, and otherwise impressing people at Mathworks enough for them to integrate MCT into Matlab.
latest version for Linux
Is the version 4.9.1.14792 available for linux?
I am able to download only version 4.9.0.14753
execution time not proportional to array length
Pavel,
Doubling array size quadruples execution time.
I'm not asking for help debugging code; this is something you need to see.
Look at it under the Matlab Profiler.
Jon
______________________________________________________________________
%MP Toolbox timing anomaly: twice the array size takes four times the execution time
clearvars
Fs = mp('48000');
freq = mp('997');
NumSecs = mp('1'); %record duration in seconds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = NumSecs*Fs;
disp(['generating ' num2str(double(NumSecs)) '-second D/A signal'])
y1 = zeros(N,1,'mp');
y2 = zeros(N,1,'mp');
cosom = cos(mp('2*pi')*freq/Fs);
tic
for n=2:N %for loop necessary in production code
y1(n) = cosom*y1(n-1); %7.5 seconds in Matlab Profiler, AMD 3995WX CPU Windows 10
y2(n) = cosom*y2(n-1); %7.5 seconds in Matlab Profiler, AMD 3995WX CPU Windows 10
end
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% DOUBLING number of array elements takes FOUR times the previous execution time %%%
NumSecs = mp('2'); %record duration in seconds
%%%%%%%%%%%%%%%%% EXACT DUPLICATE (cut & paste) of code block above %%%%%%%%%%%%%%%%
N = NumSecs*Fs;
disp(['generating ' num2str(double(NumSecs)) '-second D/A signal'])
y1 = zeros(N,1,'mp');
y2 = zeros(N,1,'mp');
cosom = cos(mp('2*pi')*freq/Fs);
tic
for n=2:N %for loop necessary in production code
y1(n) = cosom*y1(n-1); %37 seconds in Matlab Profiler, AMD 3995WX CPU Windows 10
y2(n) = cosom*y2(n-1); %37 seconds in Matlab Profiler, AMD 3995WX CPU Windows 10
end
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subroutine precision
Quoting this page:
https://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/"If input parameters a and b are multiprecision type, then all subsequent calculations inside the function will be conducted with high precision,..."
Re: Matlab 2021a under Windows 10:
User-written subroutines are executed in quadruple precision (34 digits) regardless of number of digits possessed by arguments passed to the subroutine.
For example: If arg1 has 100-digit precision and is passed to subroutine1(arg1), then computations within subroutine1 will be performed in 34-digit precision unless the following statement appears inside:
mp.Digits(mp.Digits(arg1));
R2022a incompatibility
>> mp.Test
nthroot() :Error using char/eps
Too many input arguments.
Error in nthroot (line 35)
m = x ~= 0 & (abs(x) < 1./eps('like',y)) & isfinite(n);
Error in mptest (line 1298)
d = nthroot(X,Y)-nthroot(mp(X),mp(Y));
Error in mp.Test (line 300)
mptest(); % Run tests
I see, this was helpful, thank you!
Now situation is clear. MATLAB R2022a introduced incompatibility with previous versions of MATLAB.
The issue is - R2022a chooses wrong function overload. MCT provides custom overload for 'eps' function which accepts only one argument. Two-argument calls must be re-directed to the built-in 'eps'. All previous versions of MATLAB handled this correctly. But starting from R2022a, MATLAB (incorrectly) tries to use one-argument overload for two-argument calls.
Obviously this leads to error "Too many input arguments".
Actually this issue is a bug. Overloads look-up rules are quite clear and strict in the world of OOP, it cannot be changed on a whim. Lots of existing MATLAB code might be broken because of this change.
We will add workaround for this issue to MCT, so that our overload will handle one and two-argument calls.
At the moment, please replace the content of [toolbox_dir]/@char/eps.m with the following code:
function r = eps( varargin )
if strcmpi('mp',varargin{end}) || strcmpi('mp',class(varargin{end}))
r = mpimpl(2000);
else
r = builtin('eps',varargin{:});
end
end
Let me know if this fixes the issue.
mp.read mp.write multidimensional matrix
Matlab assumes matrices to be multi-dimensional.
But Multiprecision Toolbox assumes exported matrices are only two-dimensional (can be written on a piece of paper); e.g,
Matlab input:
A = rand(3,3,4,'mp');
mp.write(A, 'mpmatrix.txt');
B = mp.read('mpmatrix.txt');
Matlab output:
Error using mp (line 1331)
Error: uneven number of elements in the rows.
Error in mp.read (line 1017)
A = mp(['[',s,']']);
This cripples communication by file transfer with Mathematica which understands Matlab .mat file format.
@"Change .txt to .mat in above; i.e"
This changes only the extension in file name (txt to mat), but file is still textual.
That is why Mathematica can read it.
Functions mp.write and mp.read were provided solely for the purpose of transferring toolbox's arrays to other software. This can be efficiently done only in textual format. In txt format we store only 2D slices for maximum compatibility with other programs.
The MATLAB's MAT files are of binary format, it is different beast. Toolbox's multiprecision arrays can be stored in binary MAT-files and loaded by toolbox on different OS - using MATLAB's standard commands, save/load.
But Mathematica cannot read it (because there is no support for our toolbox there).
Poor MCT performance on simple arithmetic
I have the following piece of code which takes 100 times longer than simple double precision calculation on my laptop. Is this expected?
s = 0.00001;
delta = 0.9;
e = linspace(0,2,1015);
gamma = 1e-05;
mp.Digits(50);
tic;
for i = 1:100
m = mp((delta/s)^2 - 1); b = mp(-2*(1i*e-gamma)./s);
p1 = m - b.^2/8; p2 = - b*(m + 2)/2;
h1 = (b.^2/4 + m).^2;
h2 = b.^6/32+(3*b.^4*m)/8+(3*b.^2*m^2)/2+27*b.^2*m+27*b.^2+2*m^3;
S = (1/2*(h2+sqrt(h2.^2-4*h1.^3))).^(1/3);
S = 1/2*sqrt(-2/3*p1+1/3*(S+h1./S));
sol = -b./4+S+1/2*sqrt(-4*S.^2-2*p1-p2./S);
sol = double(sol);
end
toc;
P.S. this is the application of Cardano formula to find the roots of a quartic polynomial when coefficients vary.
I have bolded the parts that need to be changed if run using MATLAB double precision arithmetic.
By the way, you need to compute constants (1/3, 2/3) using extended precision: mp('1/3'), mp('2/3').
Otherwise MATLAB evaluates them using double precision. Of course, this must be pre-calculated once before the loop (same as m).
There is still a room for improvement in your code.
amd threadripper
Regarding Advanpix Multiprecision Toolbox, in particular, would there be any disadvantage at all to its operation on an AMD cpu?
Customer support service by UserEcho
