0
Fixed

spdiags with low number of digits

michaloutrata 3 months ago updated by Pavel Holoborodko 2 months ago 8

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.

+1
Under review

Thank you for the report!


It is an interesting and tricky bug. Usual spdiags uses double arrays to store both - the indices and the numeric values.

If we convert the code to support mp-type, then mp array is used to store the indices and numeric values too.


Everything is Ok as long as precision is capable of storing the indices accurately. 

In the case of mp.Digits(1) this is no longer valid - and indices becomes inaccurate after hitting two digits.


I have sent you the fixed mpspdiags code via email.


HI Pavel, thanks for the very quick reaction and reply!

The issue seems to be sorted out but a different one appeared - now seemingly independet of the size-and-mp.Digits interaction. Running

mp.Digits(1);
x = sparse(1:10,ones(10,1),rand(10,1,'mp')); %%% here, the error also appears for a simple x = rand(10,'mp'); but looks a bit different and since I'm primarily working with sparse matrices, I used the relevant setup
mtrx = spdiags(ones(10,1,'mp'),0,x);

produces the error

"Error using mp/cat
Dimensions of matrices being concatenated are not consistent.

Error in mp/vertcat (line 1479)
[varargout{1:nargout}] = cat(1,varargin{:});

Error in mpspdiags

Error in mp/spdiags (line 4541)
[varargout{1:nargout}] = mpspdiags(varargin{:});"

+1

I have fixed the issue (see email). It was related to the same mix of indices/values in spdiags code (in different place).

The new spdiags will be included in the next toolbox update. Please use the file I sent you until then.

Overall, it is pretty common for M-language to mix the integers and floating point values in the same array.

Please let me know if you will find similar situations. 


The low precision of 1 decimal digit is quite an extreme case, you are the first user we know who uses it.

Could you share what problem are you solving?

Perfect - now everything seems to run smoothly! As for the application - I'm looking at iterative methods where I want to treat certain parts of the code in various precisions. Commonly, people are interested in the supported precisions of half, single and double but we were interested in treating the precision as just another parameter, which choice can be optimize. Hence the 1 or 2 digits are of interest as comparisons in some numerical tests supporting the theory rather than being of itnerest on their own (e.g., because we would really like to calculate certain things to only one decimal digits). In other words, (sadly) we don't have any particular application in which just one digit is the way to go.


Originally I even wanted to just skip the one digit but since the problem occured also with 2 and 3 digits (and presumable for the following ones as well), i thought it worthwhile to report it and then since you were already fixing it, it felt like poor manners not to report even the edge case :).

Fixed

Thank you very much for the bug report and for the details on your work, very interesting!

Let me know if you need my help with anything else. I am marking the bug as fixed.

>>mp.Digits(1)

returns 34.

mp.Digits with argument (.) is a command, not a query.

Why the latency on a command?

mp.Digits(N) returns current precision (before changing it to N) in order to be used in a constructs like this:

prevN=mp.Digits(newN)
% Compute with precision "newN"
...
mp.Digits(prevN) % restore the precision 


Our documentation explains it as following:

"Setup the default precision. All subsequent multiprecision calculations will be conducted using the newly specified number of decimal digits. The function returns the previously used default precision. If called without argument, the function simply returns the current default precision in use."

Let me know if it is too vague or not accurate - will try to re-word it in some way.