/*
* Copyright (c) 2002-2006 Samit Basu
*
* 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 2 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, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include <math.h>
#include <QtCore>
#include "Core.hpp"
#include "Exception.hpp"
#include "Array.hpp"
#include "Math.hpp"
#include "Malloc.hpp"
#include "SingularValueDecompose.hpp"
#include "QRDecompose.hpp"
#include "System.hpp"
#include "Sparse.hpp"
#include "LUDecompose.hpp"
#include "Class.hpp"
#include "Print.hpp"
#include "MemPtr.hpp"
#include "Parser.hpp"
#include <algorithm>
#undef max
#undef min
//!
//@Module DISP Display a Variable or Expression
//@@Section IO
//@@Usage
//Displays the result of a set of expressions. The @|disp| function
//takes a variable number of arguments, each of which is an expression
//to output:
//@[
// disp(expr1,expr2,...,exprn)
//@]
//This is functionally equivalent to evaluating each of the expressions
//without a semicolon after each.
//@@Example
//Here are some simple examples of using @|disp|.
//@<
//a = 32;
//b = 1:4;
//disp(a,b,pi)
//@>
//!
ArrayVector DispFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
int length;
Array C;
ArrayVector retval;
length = arg.size();
for (int i=0;i<length;i++) {
C = arg[i];
PrintArrayClassic(C,eval->getPrintLimit(),eval);
}
// retval.push_back(C);
return retval;
}
//!
//@Module SPARSE Construct a Sparse Matrix
//@@Section SPARSE
//@@Usage
//Creates a sparse matrix using one of several formats. The
//first creates a sparse matrix from a full matrix
//@[
// y = sparse(x).
//@]
//The second form creates a sparse matrix containing all zeros
//that is of the specified size (the sparse equivalent of
//@|zeros|).
//@[
// y = sparse(m,n)
//@]
//where @|m| and @|n| are integers. Just like the @|zeros| function,
//the sparse matrix returned is of type @|float|. The third form
//constructs a sparse matrix from the IJV syntax. It has two forms. The
//first version autosizes the sparse matrix
//@[
// y = sparse(i,j,v)
//@]
//while the second version uses an explicit size specification
//@[
// y = sparse(i,j,v,m,n)
//@]
//@@Tests
//@{ test_sparse1.m
//function x = test_sparse1
//a = [0,0,2,0;0,0,1,1;1,0,0,0];
//A = sparse(a);
//b = float(a);
//B = sparse(b);
//c = double(a);
//C = sparse(c);
//d = complex(a + i*a);
//D = sparse(d);
//e = dcomplex(d);
//E = sparse(e);
//f = complex(a);
//F = complex(A);
//g = dcomplex(b);
//G = dcomplex(B);
//h = float(d);
//H = float(D);
//x = testeq(a,A) & testeq(b,B) & testeq(c,C) & testeq(d,D) & testeq(e,E) & testeq(f,F) & testeq(g,G) & testeq(h,H);
//@}
//@{ test_sparse2.m
//function x = test_sparse2
//a = [0,0,2,0;0,0,1,1;1,0,0,0];
//b = [4,5,0,3];
//c = [0,0];
//d = [2,3];
//A = sparse(a);
//B = sparse(b);
//C = sparse(c);
//D = sparse(d);
//j = [a;c,d;b];
//J = [A;C,D;B];
//x = testeq(j,J);
//@}
//@{ test_sparse3.m
//function x = test_sparse3
//a = float([0,0,2,0;0,0,1,1;1,0,0,0]);
//b = float([4,5,0,3]);
//c = [0,0];
//d = [2,3];
//A = sparse(a);
//B = sparse(b);
//C = sparse(c);
//D = sparse(d);
//j = [a;c,d;b];
//J = [A;C,D;B];
//x = testeq(j,J);
//@}
//@{ test_sparse4.m
//function x = test_sparse4
//a = double([0,0,2,0;0,0,1,1;1,0,0,0]);
//b = [4,5,0,3];
//c = double([0,0]);
//d = [2,3];
//A = sparse(a);
//B = sparse(b);
//C = sparse(c);
//D = sparse(d);
//j = [a;c,d;b];
//J = [A;C,D;B];
//x = testeq(j,J);
//@}
//@{ test_sparse5.m
//function x = test_sparse5
//a = complex([0,0,2,0;0,0,1,1;1,0,0,0] + i*[3,0,3,0;0,0,1,0;0,3,5,2]);
//b = float([4,5,0,3]);
//c = float([0,0]);
//d = float([2+5*i,3]);
//A = sparse(a);
//B = sparse(b);
//C = sparse(c);
//D = sparse(d);
//j = [a;c,d;b];
//J = [A;C,D;B];
//x = testeq(j,J);
//@}
//@{ test_sparse6.m
//function x = test_sparse6
//a = dcomplex([0,0,2,0;0,0,1,1;1,0,0,0] + i*[3,0,3,0;0,0,1,0;0,3,5,2]);
//b = float([4,5,0,3]);
//c = double([0,0]);
//d = float([2+5*i,3]);
//A = sparse(a);
//B = sparse(b);
//C = sparse(c);
//D = sparse(d);
//j = [a;c,d;b];
//J = [A;C,D;B];
//x = testeq(j,J);
//@}
//@{ test_sparse7.m
//function x = test_sparse7
//a = [0,0,2,0,0;0,1,0,1,0;0,0,0,1,3];
//A = sparse(a);
//b = [1,4,5;2,5,15;3,6,11];
//C = A(b);
//c = a(b);
//x = testeq(c,C);
//@}
//@{ test_sparse8.m
//function x = test_sparse8
//ar = [0,0,2,0,0;0,1,0,1,0;0,0,0,1,0];
//ai = [1,0,3,0,0;0,2,3,-1,3;0,1,0,4,3];
//a = complex(ar+i*ai);
//A = sparse(a);
//b = [1,4,5;2,5,15;3,6,11];
//C = A(b);
//c = a(b);
//x = testeq(c,C);
//@}
//@{ test_sparse9.m
//function x = test_sparse9
//ar = [0,0,2,0,0;0,1,0,1,0;0,0,0,1,0];
//ai = [1,0,3,0,0;0,2,3,-1,3;0,1,0,4,3];
//a = dcomplex(ar+i*ai);
//A = sparse(a);
//b = [1,4,5;2,5,15;3,6,11];
//C = A(b);
//c = a(b);
//x = testeq(c,C);
//@}
//@{ test_sparse10.m
//function x = test_sparse10
//a = int32(10*rand(8));
//a(a<7) = 0;
//A = sparse(a);
//b = [5;2;3;3;1];
//c = [1,4,6,8];
//A = sparse(a);
//C = A(b,c);
//c = a(b,c);
//x = testeq(c,C);
//@}
//@{ test_sparse11.m
//function x = test_sparse11
//a = float(10*rand(8));
//a(a<7) = 0;
//A = sparse(a);
//b = [5;2;3;3;1];
//c = [1,4,6,8];
//A = sparse(a);
//C = A(b,c);
//c = a(b,c);
//x = testeq(c,C);
//@}
//@{ test_sparse12.m
//function x = test_sparse12
//a = double(10*rand(8));
//a(a<7) = 0;
//A = sparse(a);
//b = [5;2;3;3;1];
//c = [1,4,6,8];
//A = sparse(a);
//C = A(b,c);
//c = a(b,c);
//x = testeq(c,C);
//@}
//@{ test_sparse13.m
//function x = test_sparse13
//ai = int32(10*rand(8));
//ar = int32(10*rand(8));
//ai(ai<7) = 0;
//ar(ar<7) = 0;
//a = complex(ar+i*ai);
//A = sparse(a);
//b = [5;2;3;3;1];
//d = [1,4,6,8];
//A = sparse(a);
//C = A(b,d);
//c = a(b,d);
//x = testeq(c,C);
//@}
//@{ test_sparse14.m
//function x = test_sparse14
//ai = int32(10*rand(8));
//ar = int32(10*rand(8));
//ai(ai<7) = 0;
//ar(ar<7) = 0;
//a = dcomplex(ar+i*ai);
//A = sparse(a);
//b = [5;2;3;3;1];
//c = [1,4,6,8];
//A = sparse(a);
//C = A(b,c);
//c = a(b,c);
//x = testeq(c,C);
//@}
//@{ test_sparse15.m
//function x = test_sparse15
//a = [1,2,0,0,4;3,2,0,0,5;0,0,3,0,2];
//A = sparse(a);
//c = a(2,5);
//C = A(2,5);
//x = testeq(c,C);
//@}
//@{ test_sparse16.m
//function x = test_sparse16
//a = float([1,2,0,0,4;3,2,0,0,5;0,0,3,0,2]);
//A = sparse(a);
//c = a(2,5);
//C = A(2,5);
//x = testeq(c,C);
//@}
//@{ test_sparse17.m
//function x = test_sparse17
//a = double([1,2,0,0,4;3,2,0,0,5;0,0,3,0,2]);
//A = sparse(a);
//c = a(2,5);
//C = A(2,5);
//x = testeq(c,C);
//@}
//@{ test_sparse18.m
//function x = test_sparse18
//ai = [0,2,6,0,1;3,0,3,0,2;0,0,3,0,2];
//ar = [1,2,0,0,4;3,2,0,0,5;0,0,3,0,2];
//a = complex(ar+i*ai);
//A = sparse(a);
//c = a(2,5);
//C = A(2,5);
//x = testeq(c,C);
//@}
//@{ test_sparse19.m
//function x = test_sparse19
//ai = [0,2,6,0,1;3,0,3,0,2;0,0,3,0,2];
//ar = [1,2,0,0,4;3,2,0,0,5;0,0,3,0,2];
//a = dcomplex(ar+i*ai);
//A = sparse(a);
//c = a(2,5);
//C = A(2,5);
//x = testeq(c,C);
//@}
//@{ test_sparse22.m
//function x = test_sparse22
//a = randn(1000);
//A = sparse(a);
//rc = int32(999*rand(7000,1))+1;
//cc = int32(999*rand(7000,1))+1;
//lc = rc+1000*(cc-1);
//b = a(lc);
//B = A(lc);
//x = testeq(b,B);
//@}
//@{ test_sparse23.m
//function x = test_sparse23
//N = 1000;
//M = 7000;
//a = double(40*randn(N));
//A = sparse(a);
//rc = int32((N-1)*rand(M,1))+1;
//cc = int32((N-1)*rand(M,1))+1;
//lc = rc+N*(cc-1);
//b = a(lc);
//B = A(lc);
//x = testeq(b,B);
//@}
//@{ test_sparse24.m
//function x = test_sparse24
//N = 10;
//M = 70;
//a = double(40*randn(N));
//A = sparse(a);
//rc = int32((N-1)*rand(M,1))+1;
//cc = int32((N-1)*rand(M,1))+1;
//lc = rc+N*(cc-1);
//b = a(lc);
//B = A(lc);
//x = testeq(b,B);
//@}
//@{ test_sparse25.m
//function x = test_sparse25
//N = 1000;
//M = 7000;
//a = complex(40*randn(N)+32*randn(N));
//A = sparse(a);
//rc = int32((N-1)*rand(M,1))+1;
//cc = int32((N-1)*rand(M,1))+1;
//lc = rc+N*(cc-1);
//b = a(lc);
//B = A(lc);
//x = testeq(b,B);
//@}
//@{ test_sparse26.m
//function x = test_sparse26
//N = 1000;
//M = 7000;
//a = dcomplex(40*randn(N)+32*randn(N));
//A = sparse(a);
//rc = int32((N-1)*rand(M,1))+1;
//cc = int32((N-1)*rand(M,1))+1;
//lc = rc+N*(cc-1);
//b = a(lc);
//B = A(lc);
//x = testeq(b,B);
//@}
//@{ test_sparse27.m
//function x = test_sparse27
//a = [1,0,3,4,5;6,2,3,5,0;0,0,0,0,2];
//A = sparse(a);
//p = [3;4;5;9;10];
//a(p) = 7;
//A(p) = 7;
//x = testeq(a,A);
//@}
//@{ test_sparse28.m
//function x = test_sparse28
//a = [1,0,3,4,5;6,2,3,5,0;0,0,0,0,2];
//A = sparse(a);
//p = [3;4;5;9;10];
//a(p) = [4,6,8,0,3];
//A(p) = [4,6,8,0,3];
//x = testeq(a,A);
//@}
//@{ test_sparse29.m
//function x = test_sparse29
//a = [1,0,3,4,5;6,2,3,5,0;0,0,0,0,2];
//A = sparse(a);
//p = [3;4;5;9;10;20];
//a(p) = [4,6,8,0,3,7];
//A(p) = [4,6,8,0,3,7];
//x = testeq(a',A);
//@}
//@{ test_sparse30.m
//function x = test_sparse30
//ar = [1,0,3,4,5;6,2,3,5,0;0,0,0,0,2];
//ai = [0,0,2,0,4;1,0,2,3,0;1,0,0,0,3];
//a = complex(ar+ai*i);
//A = sparse(a);
//p = [3;4;5;9;10];
//a(p) = complex(7+2*i);
//A(p) = complex(7+2*i);
//x = testeq(a,A);
//@}
//@{ test_sparse31.m
//function x = test_sparse31
//a = [1,0,3,4,5;6,2,3,5,0;0,0,0,0,2;0,4,5,0,3];
//A = sparse(a);
//b = [3;2];
//c = [2,4,3];
//d = [9,3,2,4,9,0];
//a(b,c) = d;
//A(b,c) = d;
//x = testeq(a,A);
//@}
//@{ test_sparse32.m
//function x = test_sparse32
//a = [1,0,3,4,5;6,2,3,5,0;0,0,0,0,2;0,4,5,0,3];
//A = sparse(a);
//b = [3;2;5];
//c = [2,4,3];
//d = [9,3,2,4,9,0,5,4,2];
//a(b,c) = d;
//A(b,c) = d;
//x = testeq(a,A);
//@}
//@{ test_sparse33.m
//function x = test_sparse33
//ar = [1,0,3,4,5;6,2,3,5,0;0,0,0,0,2;0,4,5,0,3];
//ai = [0,2,9,0,2;0,0,0,2,0;1,3,4,0,2;1,0,0,0,2];
//a = complex(ar+ai);
//A = sparse(a);
//b = [3;2;5];
//c = [2,4,3];
//d = complex([9,3,2,4,9,0,5,4,2] + 3*i);
//a(b,c) = d;
//A(b,c) = d;
//x = testeq(a,A);
//@}
//@{ test_sparse34.m
//function x = test_sparse34
//a = sparse([0;0;2;3;0;0;4;5;0;0;0;2]);
//n = [7,9,1,3,2,3,1,2,2];
//b = full(a);
//a(n) = [];
//b(n) = [];
//x = testeq(a,b);
//@}
//@{ test_sparse35.m
//function x = test_sparse35
//s = [0,3,4,5;0,0,0,0;0,2,3,0];
//a = sparse(s);
//x = testeq(full(a),s) & testeq(full(float(a)),s) & testeq(full(double(a)),s) & testeq(full(complex(a)),s) & testeq(full(dcomplex(a)),s);
//@}
//@{ test_sparse36.m
//function x = test_sparse36
//a = rand(8) + i*rand(8);
//A = sparse(a);
//B = float(A);
//b = full(B);
//x = testeq(b,B);
//@}
//@{ sparse_test_mat.m
//function [y,z] = sparse_test_mat(typ, rows, cols)
// if (~exist('cols')) cols = rows; end
// a = rand(rows,cols);
// a(a>0.1) = 0;
// b = rand(rows,cols);
// b(b>0.1) = 0;
// switch(typ)
// case 'int32'
// z = (int32(100*a));
// case 'float'
// z = (float(a));
// case 'double'
// z = (double(a));
// case 'complex'
// z = (float(a) + i*float(b));
// case 'dcomplex'
// z = (a + i*b);
// end
// y = sparse(z);
//@}
//@{ test_sparse62.m
//% Test sparse matrix array constructor
//function x = test_sparse62
//[yi1,zi1] = sparse_test_mat('int32',300,400);
//[yf1,zf1] = sparse_test_mat('float',300,400);
//[yd1,zd1] = sparse_test_mat('double',300,400);
//[yc1,zc1] = sparse_test_mat('complex',300,400);
//[yz1,zz1] = sparse_test_mat('dcomplex',300,400);
//[yi2,zi2] = sparse_test_mat('int32',300,450);
//[yf2,zf2] = sparse_test_mat('float',300,450);
//[yd2,zd2] = sparse_test_mat('double',300,450);
//[yc2,zc2] = sparse_test_mat('complex',300,450);
//[yz2,zz2] = sparse_test_mat('dcomplex',300,450);
//[yi3,zi3] = sparse_test_mat('int32',30,400);
//[yf3,zf3] = sparse_test_mat('float',30,400);
//[yd3,zd3] = sparse_test_mat('double',30,400);
//[yc3,zc3] = sparse_test_mat('complex',30,400);
//[yz3,zz3] = sparse_test_mat('dcomplex',30,400);
//[yi4,zi4] = sparse_test_mat('int32',30,450);
//[yf4,zf4] = sparse_test_mat('float',30,450);
//[yd4,zd4] = sparse_test_mat('double',30,450);
//[yc4,zc4] = sparse_test_mat('complex',30,450);
//[yz4,zz4] = sparse_test_mat('dcomplex',30,450);
//a1 = [yi1,yi2;yi3,yi4];
//b1 = [zi1,zi2;zi3,zi4];
//a2 = [yf1,yf2;yf3,yf4];
//b2 = [zf1,zf2;zf3,zf4];
//a3 = [yd1,yd2;yd3,yd4];
//b3 = [zd1,zd2;zd3,zd4];
//a4 = [yc1,yc2;yc3,yc4];
//b4 = [zc1,zc2;zc3,zc4];
//a5 = [yz1,yz2;yz3,yz4];
//b5 = [zz1,zz2;zz3,zz4];
//x = testeq(a1,b1) & testeq(a2,b2) & testeq(a3,b3) & testeq(a4,b4) & testeq(a5,b5);
//@}
//!
ArrayVector SparseFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
if (arg.size() == 1) {
Array r(arg[0]);
if ((r.dataClass() != FM_LOGICAL) && (r.dataClass() < FM_INT32))
r.promoteType(FM_INT32);
r.makeSparse();
return singleArrayVector(r);
} else if (arg.size() == 2) {
Array m_arg(arg[0]);
Array n_arg(arg[1]);
int m, n;
m = m_arg.getContentsAsIntegerScalar();
n = n_arg.getContentsAsIntegerScalar();
Dimensions dim(m,n);
return singleArrayVector(Array(FM_FLOAT,dim,SparseFloatZeros(m,n),true));
} else if (arg.size() == 3) {
Array i_arg(arg[0]);
Array j_arg(arg[1]);
Array v_arg(arg[2]);
i_arg.promoteType(FM_UINT32);
j_arg.promoteType(FM_UINT32);
if (v_arg.dataClass() != FM_LOGICAL && v_arg.dataClass() < FM_INT32)
v_arg.promoteType(FM_INT32);
int ilen, jlen, vlen;
ilen = i_arg.getLength();
jlen = j_arg.getLength();
vlen = v_arg.getLength();
int olen;
olen = ilen > jlen ? ilen : jlen;
olen = vlen > olen ? vlen : olen;
int istride, jstride, vstride;
if (olen > 1) {
if (ilen == 1)
istride = 0;
else if (ilen == olen)
istride = 1;
else
throw Exception("in I, J, V format, all three vectors must be the same size or be scalars");
if (jlen == 1)
jstride = 0;
else if (jlen == olen)
jstride = 1;
else
throw Exception("in I, J, V format, all three vectors must be the same size or be scalars");
if (vlen == 1)
vstride = 0;
else if (vlen == olen)
vstride = 1;
else
throw Exception("in I, J, V format, all three vectors must be the same size or be scalars");
}
// Calculate the number of rows in the matrix
uint32 *ip = (uint32*) i_arg.getDataPointer();
int rows = 0;
for (int i=0;i<ilen;i++)
rows = (ip[i] > rows) ? ip[i] : rows;
uint32 *jp = (uint32*) j_arg.getDataPointer();
int cols = 0;
for (int j=0;j<jlen;j++)
cols = (jp[j] > cols) ? jp[j] : cols;
Dimensions dim(rows,cols);
return singleArrayVector(Array(v_arg.dataClass(),dim,
makeSparseFromIJV(v_arg.dataClass(),
rows,cols,olen,
ip,istride,jp,jstride,
v_arg.getDataPointer(),
vstride),
true));
} else if (arg.size() >= 5) {
if (arg.size() > 5)
eval->warningMessage("extra arguments to sparse (nnz to reserve) ignored");
Array i_arg(arg[0]);
Array j_arg(arg[1]);
Array v_arg(arg[2]);
i_arg.promoteType(FM_UINT32);
j_arg.promoteType(FM_UINT32);
if (v_arg.dataClass() != FM_LOGICAL && v_arg.dataClass() < FM_INT32)
v_arg.promoteType(FM_INT32);
int ilen, jlen, vlen;
ilen = i_arg.getLength();
jlen = j_arg.getLength();
vlen = v_arg.getLength();
int olen;
olen = ilen > jlen ? ilen : jlen;
olen = vlen > olen ? vlen : olen;
int istride, jstride, vstride;
if (olen > 1) {
if (ilen == 1)
istride = 0;
else if (ilen == olen)
istride = 1;
else
throw Exception("in I, J, V format, all three vectors must be the same size or be scalars");
if (jlen == 1)
jstride = 0;
else if (jlen == olen)
jstride = 1;
else
throw Exception("in I, J, V format, all three vectors must be the same size or be scalars");
if (vlen == 1)
vstride = 0;
else if (vlen == olen)
vstride = 1;
else
throw Exception("in I, J, V format, all three vectors must be the same size or be scalars");
}
Array m_arg(arg[3]);
Array n_arg(arg[4]);
int rows = m_arg.getContentsAsIntegerScalar();
int cols = n_arg.getContentsAsIntegerScalar();
Dimensions dim(rows,cols);
uint32 *ip = (uint32*) i_arg.getDataPointer();
uint32 *jp = (uint32*) j_arg.getDataPointer();
return singleArrayVector(Array(v_arg.dataClass(),dim,
makeSparseFromIJV(v_arg.dataClass(),
rows,cols,olen,
ip,istride,jp,jstride,
v_arg.getDataPointer(),
vstride),
true));
}
throw Exception("unrecognized form of sparse - see help for the allowed forms of sparse");
}
//!
//@Module INV Invert Matrix
//@@Section TRANSFORMS
//@@Usage
//Inverts the argument matrix, provided it is square and invertible.
//The syntax for its use is
//@[
// y = inv(x)
//@]
//Internally, the @|inv| function uses the matrix divide operators.
//For sparse matrices, a sparse matrix solver is used.
//@@Example
//Here we invert some simple matrices
//@<
//a = randi(zeros(3),5*ones(3))
//b = inv(a)
//a*b
//b*a
//@>
//!
ArrayVector InvFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
if (arg.size() != 1)
throw Exception("inv function needs at least one argument");
Array r(arg[0]);
return singleArrayVector(InvertMatrix(r,eval));
}
//!
//@Module FULL Convert Sparse Matrix to Full Matrix
//@@Section SPARSE
//@@Usage
//Converts a sparse matrix to a full matrix. The syntax for
//its use is
//@[
// y = full(x)
//@]
//The type of @|x| is preserved. Be careful with the function.
//As a general rule of thumb, if you can work with the @|full|
//representation of a function, you probably do not need the
//sparse representation.
//@@Example
//Here we convert a full matrix to a sparse one, and back again.
//@<
//a = [1,0,4,2,0;0,0,0,0,0;0,1,0,0,2]
//A = sparse(a)
//full(A)
//@>
//@@Tests
//@{ test_sparse57.m
//% Test makeDenseArray function
//function x = test_sparse57
//[yi,zi] = sparse_test_mat('int32',400);
//[yf,zf] = sparse_test_mat('float',400);
//[yd,zd] = sparse_test_mat('double',400);
//[yc,zc] = sparse_test_mat('complex',400);
//[yz,zz] = sparse_test_mat('dcomplex',400);
//x = testeq(full(yi),zi) & testeq(full(yf),zf) & testeq(full(yd),zd) & testeq(full(yc),zc) & testeq(full(yz),zz);
//@}
//!
ArrayVector FullFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 1)
throw Exception("Need one argument to full function");
Array r(arg[0]);
r.makeDense();
return singleArrayVector(r);
}
//!
//@Module LU LU Decomposition for Matrices
//@@Section TRANSFORMS
//@@Usage
//Computes the LU decomposition for a matrix. The form of the
//command depends on the type of the argument. For full (non-sparse)
//matrices, the primary form for @|lu| is
//@[
// [L,U,P] = lu(A),
//@]
//where @|L| is lower triangular, @|U| is upper triangular, and
//@|P| is a permutation matrix such that @|L*U = P*A|. The second form is
//@[
// [V,U] = lu(A),
//@]
//where @|V| is @|P'*L| (a row-permuted lower triangular matrix),
//and @|U| is upper triangular. For sparse, square matrices,
//the LU decomposition has the following form:
//@[
// [L,U,P,Q,R] = lu(A),
//@]
//where @|A| is a sparse matrix of either @|double| or @|dcomplex| type.
//The matrices are such that @|L*U=P*R*A*Q|, where @|L| is a lower triangular
//matrix, @|U| is upper triangular, @|P| and @|Q| are permutation vectors
//and @|R| is a diagonal matrix of row scaling factors. The decomposition
// is computed using UMFPACK for sparse matrices, and LAPACK for dense
// matrices.
//@@Example
//First, we compute the LU decomposition of a dense matrix.
//@<
//a = float([1,2,3;4,5,8;10,12,3])
//[l,u,p] = lu(a)
//l*u
//p*a
//@>
//Now we repeat the exercise with a sparse matrix, and demonstrate
//the use of the permutation vectors.
//@<
//a = sparse([1,0,0,4;3,2,0,0;0,0,0,1;4,3,2,4])
//[l,u,p,q,r] = lu(a)
//full(l*a)
//b = r*a
//full(b(p,q))
//@>
//@@Tests
//@{ test_lu1.m
//% Test the LU decomposition with 2 argument return for full matrices
//function x = test_lu1
//A = float(randn(20));
//[L,U] = lu(A);
//x = testeq(100,L*U,A);
//A = double(randn(20));
//[L,U] = lu(A);
//x = x & testeq(100,L*U,A);
//A = complex(randn(20)+i*randn(20));
//[L,U] = lu(A);
//x = x & testeq(100,L*U,A);
//A = dcomplex(randn(20)+i*randn(20));
//[L,U] = lu(A);
//x = x & testeq(100,L*U,A);
//
//function x = testeq(tmag,a,b)
// d = full(a)-full(b);
// if (strcmp(typeof(d),'double') | strcmp(typeof(d),'dcomplex'))
// x = isempty(find(abs(d)>tmag*eps));
// else
// x = isempty(find(abs(d)>tmag*feps));
// end
//@}
//@{ test_lu2.m
//% Test the LU decomposition with 3 argument return for full matrices
//function x = test_lu2
//A = float(randn(20));
//[L,U,P] = lu(A);
//x = testeq(100,L*U,P*A);
//A = double(randn(20));
//[L,U,P] = lu(A);
//x = x & testeq(100,L*U,P*A);
//A = complex(randn(20)+i*randn(20));
//[L,U,P] = lu(A);
//x = x & testeq(100,L*U,P*A);
//A = dcomplex(randn(20)+i*randn(20));
//[L,U,P] = lu(A);
//x = x & testeq(100,L*U,P*A);
//
//function x = testeq(tmag,a,b)
// d = full(a)-full(b);
// if (strcmp(typeof(d),'double') | strcmp(typeof(d),'dcomplex'))
// x = isempty(find(abs(d)>tmag*eps));
// else
// x = isempty(find(abs(d)>tmag*feps));
// end
//@}
//@{ test_sparse75.m
//% Test sparse matrix array finite test
//function x = test_sparse75
//[yi1,zi1] = sparse_test_mat('int32',300,400);
//[yf1,zf1] = sparse_test_mat('float',300,400);
//[yd1,zd1] = sparse_test_mat('double',300,400);
//[yc1,zc1] = sparse_test_mat('complex',300,400);
//[yz1,zz1] = sparse_test_mat('dcomplex',300,400);
//x = 1;
//yi1(150,200) = inf;
//try
// [l,u] = lu(yi1);
// x = 0;
//catch
// x = x & 1;
//end
//try
// [l,u] = lu(yf1);
// x = 0;
//catch
// x = x & 1;
//end
//try
// [l,u] = lu(yd1);
// x = 0;
//catch
// x = x & 1;
//end
//try
// [l,u] = lu(yc1);
// x = 0;
//catch
// x = x & 1;
//end
//try
// [l,u] = lu(yz1);
// x = 0;
//catch
// x = x & 1;
//end
//@}
//!
ArrayVector LUFunction(int nargout, const ArrayVector& arg) {
if (arg.size() < 1)
throw Exception("lu function requires at least one argument - the matrix to decompose.");
Array A(arg[0]);
if (A.isReferenceType())
throw Exception("cannot apply lu decomposition to reference types.");
if (!A.is2D())
throw Exception("cannot apply matrix operations to N-dimensional arrays");
if (!A.sparse()) {
if (A.anyNotFinite())
throw Exception("lu decomposition only defined for matrices with finite entries.");
return LUDecompose(nargout,A);
} else {
return SparseLUDecompose(nargout,A);
}
return ArrayVector();
}
//!
//@Module GETLINE Get a Line of Input from User
//@@Section IO
//@@Usage
//Reads a line (as a string) from the user. This function has
//two syntaxes. The first is
//@[
// a = getline(prompt)
//@]
//where @|prompt| is a prompt supplied to the user for the query.
//The second syntax omits the @|prompt| argument:
//@[
// a = getline
//@]
//Note that this function requires command line input, i.e., it
//will only operate correctly for programs or scripts written
//to run inside the FreeMat GUI environment or from the X11 terminal.
//If you build a stand-alone application and expect it to operate
//cross-platform, do not use this function (unless you include
//the FreeMat console in the final application).
//!
ArrayVector GetLineFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
string prompt;
if (arg.size() < 1)
prompt = "";
else {
Array A(arg[0]);
if (!A.isString())
throw Exception("getline requires a string prompt");
prompt = A.getContentsAsString();
}
return singleArrayVector(Array::stringConstructor(eval->getLine(prompt)));
}
ArrayVector GenEigFunction(int nargout, const ArrayVector &arg, Interpreter* m_eval) {
Array A(arg[0]);
Array B(arg[1]);
if (A.sparse() || B.sparse())
throw Exception("eig only defined for full matrices.");
if (A.anyNotFinite() || B.anyNotFinite())
throw Exception("eig only defined for matrices with finite entries.");
if (A.isReferenceType() || B.isReferenceType())
throw Exception("Cannot apply eigendecomposition to reference types.");
if (!A.is2D() || !B.is2D())
throw Exception("Cannot apply eigendecomposition to N-Dimensional arrays.");
if (A.getDimensionLength(0) != A.getDimensionLength(1))
throw Exception("Cannot eigendecompose a non-square matrix.");
if (B.getDimensionLength(0) != B.getDimensionLength(1))
throw Exception("Cannot eigendecompose a non-square matrix.");
if (A.getDimensionLength(0) != B.getDimensionLength(0))
throw Exception("B and A must be the same size when computing a generalized eigendecomposition");
// Handle the type of A - if it is an integer type, then promote to double
Class Aclass = A.dataClass();
if (Aclass < FM_FLOAT) {
A.promoteType(FM_DOUBLE);
Aclass = FM_DOUBLE;
}
Class Bclass = B.dataClass();
if (Bclass < Aclass) {
B.promoteType(Aclass);
Bclass = Aclass;
} else {
A.promoteType(Bclass);
Aclass = Bclass;
}
ArrayVector retval;
Array V, D;
if (nargout > 1) {
if (A.isSymmetric() && B.isSymmetric()) {
if (!GeneralizedEigenDecomposeFullSymmetric(A,B,V,D,m_eval))
GeneralizedEigenDecomposeFullGeneral(A,B,V,D,m_eval);
} else
GeneralizedEigenDecomposeFullGeneral(A,B,V,D,m_eval);
retval.push_back(V);
retval.push_back(D);
} else {
if (A.isSymmetric() && B.isSymmetric()) {
if (!GeneralizedEigenDecomposeCompactSymmetric(A,B,D,m_eval))
GeneralizedEigenDecomposeCompactGeneral(A,B,D,m_eval);
} else
GeneralizedEigenDecomposeCompactGeneral(A,B,D,m_eval);
retval.push_back(D);
}
return retval;
}
//!
//@Module EIG Eigendecomposition of a Matrix
//@@Section TRANSFORMS
//@@Usage
//Computes the eigendecomposition of a square matrix. The @|eig| function
//has several forms. The first returns only the eigenvalues of the matrix:
//@[
// s = eig(A)
//@]
//The second form returns both the eigenvectors and eigenvalues as two
//matrices (the eigenvalues are stored in a diagonal matrix):
//@[
// [V,D] = eig(A)
//@]
//where @|D| is the diagonal matrix of eigenvalues, and @|V| is the
//matrix of eigenvectors.
//
//Eigenvalues and eigenvectors for asymmetric matrices @|A| normally
//are computed with balancing applied. Balancing is a scaling step
//that normaly improves the quality of the eigenvalues and eigenvectors.
//In some instances (see the Function Internals section for more details)
//it is necessary to disable balancing. For these cases, two additional
//forms of @|eig| are available:
//@[
// s = eig(A,'nobalance'),
//@]
//which computes the eigenvalues of @|A| only, and does not balance
//the matrix prior to computation. Similarly,
//@[
// [V,D] = eig(A,'nobalance')
//@]
//recovers both the eigenvectors and eigenvalues of @|A| without balancing.
//Note that the 'nobalance' option has no affect on symmetric matrices.
//
//FreeMat also provides the ability to calculate generalized eigenvalues
//and eigenvectors. Similarly to the regular case, there are two forms
//for @|eig| when computing generalized eigenvector (see the Function
//Internals section for a description of what a generalized eigenvector is).
//The first returns only the generalized eigenvalues of the matrix
//pair @|A,B|
//@[
// s = eig(A,B)
//@]
//The second form also computes the generalized eigenvectors, and is
//accessible via
//@[
// [V,D] = eig(A,B)
//@]
//@@Function Internals
//Recall that @|v| is an eigenvector of @|A| with associated eigenvalue
//@|d| if
//\[
// A v = d v.
//\]
//This decomposition can be written in matrix form as
//\[
// A V = V D
//\]
//where
//\[
// V = [v_1,v_2,\ldots,v_n], D = \mathrm{diag}(d_1,d_2,\ldots,d_n).
//\]
//The @|eig| function uses the @|LAPACK| class of functions @|GEEVX|
//to compute the eigenvalue decomposition for non-symmetric
//(or non-Hermitian) matrices @|A|. For symmetric matrices, @|SSYEV|
// and @|DSYEV| are used for @|float| and @|double| matrices (respectively).
//For Hermitian matrices, @|CHEEV| and @|ZHEEV| are used for @|complex|
//and @|dcomplex| matrices.
//
//For some matrices, the process of balancing (in which the rows and
//columns of the matrix are pre-scaled to facilitate the search for
//eigenvalues) is detrimental to the quality of the final solution.
//This is particularly true if the matrix contains some elements on
//the order of round off error. See the Example section for an example.
//
//A generalized eigenvector of the matrix pair @|A,B| is simply a
//vector @|v| with associated eigenvalue @|d| such that
//\[
// A v = d B v,
//\]
//where @|B| is a square matrix of the same size as @|A|. This
//decomposition can be written in matrix form as
//\[
// A V = B V D
//\]
//where
//\[
// V = [v_1,v_2,\ldots,v_n], D = \mathrm{diag}(d_1,d_2,\ldots,d_n).
//\]
//For general matrices @|A| and @|B|, the @|GGEV| class of routines are
//used to compute the generalized eigendecomposition. If howevever,
//@|A| and @|B| are both symmetric (or Hermitian, as appropriate),
//Then FreeMat first attempts to use @|SSYGV| and @|DSYGV| for @|float|
//and @|double| arguments and @|CHEGV| and @|ZHEGV| for @|complex|
//and @|dcomplex| arguments (respectively). These routines requires
//that @|B| also be positive definite, and if it fails to be, FreeMat
//will revert to the routines used for general arguments.
//@@Example
//Some examples of eigenvalue decompositions. First, for a diagonal
//matrix, the eigenvalues are the diagonal elements of the matrix.
//@<
//A = diag([1.02f,3.04f,1.53f])
//eig(A)
//@>
//Next, we compute the eigenvalues of an upper triangular matrix,
//where the eigenvalues are again the diagonal elements.
//@<
//A = [1.0f,3.0f,4.0f;0,2.0f,6.7f;0.0f,0.0f,1.0f]
//eig(A)
//@>
//Next, we compute the complete eigenvalue decomposition of
//a random matrix, and then demonstrate the accuracy of the solution
//@<
//A = float(randn(2))
//[V,D] = eig(A)
//A*V - V*D
//@>
//Now, we consider a matrix that requires the nobalance option
//to compute the eigenvalues and eigenvectors properly. Here is
//an example from MATLAB's manual.
//@<
//B = [3,-2,-.9,2*eps;-2,4,1,-eps;-eps/4,eps/2,-1,0;-.5,-.5,.1,1]
//[VB,DB] = eig(B)
//B*VB - VB*DB
//[VN,DN] = eig(B,'nobalance')
//B*VN - VN*DN
//@>
//@@Tests
//@{ test_eig1.m
//% Test eigenvalue function - general matrices
//function t = test_eig1
//% First the float version
//t1all = 1;
//for i=2:4:100
// a = float(randn(i));
// [v,d] = eig(a);
// emat = a*v - v*d;
// er = max(abs(emat(:)));
// bnd = 100*max(diag(abs(d)))*feps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('full: compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t1all = t1all & t1 & tb;
//end
//% Now the double version
//t2all = 1;
//for i=2:4:100
// a = double(randn(i));
// [v,d] = eig(a);
// emat = a*v - v*d;
// er = max(abs(emat(:)));
// bnd = 100*max(diag(abs(d)))*eps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('double: compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t2all = t2all & t1 & tb;
//end
//% Now the complex version
//t3all = 1;
//for i=2:4:100
// a = complex(randn(i)+j*randn(i));
// [v,d] = eig(a);
// emat = a*v - v*d;
// er = max(abs(emat(:)));
// bnd = 100*max(diag(abs(d)))*feps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('complex: compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t3all = t3all & t1 & tb;
//end
//% Now the double version
//t4all = 1;
//for i=2:4:100
// a = dcomplex(randn(i)+j*randn(i));
// [v,d] = eig(a);
// emat = a*v - v*d;
// er = max(abs(emat(:)));
// bnd = 100*max(diag(abs(d)))*eps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('dcomplex: compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t4all = t4all & t1 & tb;
//end
//t = t1all & t2all & t3all & t4all;
//@}
//@{ test_eig2.m
//% Test eigenvalue function - symmetric matrices
//function t = test_eig2
//% First the float version
//t1all = 1;
//for i=2:4:100
// a = float(randn(i)); a = a + a';
// [v,d] = eig(a);
// emat = a*v - v*d;
// er = max(abs(emat(:)));
// bnd = 10*max(diag(abs(d)))*feps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t1all = t1all & t1 & tb;
//end
//% Now the double version
//t2all = 1;
//for i=2:4:100
// a = double(randn(i)); a = a + a';
// [v,d] = eig(a);
// emat = a*v - v*d;
// er = max(abs(emat(:)));
// bnd = 10*max(diag(abs(d)))*eps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t2all = t2all & t1 & tb;
//end
//% Now the complex version
//t3all = 1;
//for i=2:4:100
// a = complex(randn(i)+j*randn(i)); a = a + a';
// [v,d] = eig(a);
// emat = a*v - v*d;
// er = max(abs(emat(:)));
// bnd = 10*max(diag(abs(d)))*feps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t3all = t3all & t1 & tb;
//end
//% Now the double version
//t4all = 1;
//for i=2:4:100
// a = dcomplex(randn(i)+j*randn(i)); a = a + a';
// [v,d] = eig(a);
// emat = a*v - v*d;
// er = max(abs(emat(:)));
// bnd = 10*max(diag(abs(d)))*eps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t4all = t4all & t1 & tb;
//end
//t = t1all & t2all & t3all & t4all;
//@}
//@{ test_eig3.m
//% Test the 'nobalance' option for eig
//function t = test_eig3
//B = [3,-2,-.9,2*eps;-2,4,1,-eps;-eps/4,eps/2,-1,0;-.5,-.5,.1,1];
//[VN,DN] = eig(B,'nobalance');
//er = B*VN - VN*DN;
//er = max(abs(er(:)));
//bnd = 1.2*max(diag(DN))*eps*8;
//t = (er < bnd);
//@}
//@{ test_eig4.m
//% Test generalized eigenvalue function - general matrices
//function t = test_eig4
//% First the float version
//t1all = 1;
//for i=2:4:100
// a = float(randn(i));
// b = float(randn(i));
// [v,d] = eig(a,b);
// emat = a*v - b*v*d;
// er = max(abs(emat(:)));
// bnd = 8*max(diag(abs(d)))*feps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a,b);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('float: compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t1all = t1all & t1 & tb;
//end
//% Now the double version
//t2all = 1;
//for i=2:4:100
// a = double(randn(i));
// b = double(randn(i));
// [v,d] = eig(a,b);
// emat = a*v - b*v*d;
// er = max(abs(emat(:)));
// bnd = 8*max(diag(abs(d)))*eps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a,b);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('double: compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t2all = t2all & t1 & tb;
//end
//% Now the complex version
//t3all = 1;
//for i=2:4:100
// a = complex(randn(i)+j*randn(i));
// b = complex(randn(i)+j*randn(i));
// [v,d] = eig(a,b);
// emat = a*v - b*v*d;
// er = max(abs(emat(:)));
// bnd = 8*max(diag(abs(d)))*feps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a,b);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('complex: compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t3all = t3all & t1 & tb;
//end
//% Now the double version
//t4all = 1;
//for i=2:4:100
// a = dcomplex(randn(i)+j*randn(i));
// b = dcomplex(randn(i)+j*randn(i));
// [v,d] = eig(a,b);
// emat = a*v - b*v*d;
// er = max(abs(emat(:)));
// bnd = 8*max(diag(abs(d)))*eps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a,b);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('dcomplex: compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t4all = t4all & t1 & tb;
//end
//t = t1all & t2all & t3all & t4all;
//@}
//@{ test_eig5.m
//% Test generalized eigenvalue function - symmetric matrices
//function t = test_eig5
//% First the float version
//t1all = 1;
//for i=2:4:100
// a = float(randn(i)); a = a + a';
// b = float(randn(i)); b = b*b';
// [v,d] = eig(a,b);
// emat = a*v - b*v*d;
// er = max(abs(emat(:)));
// bnd = 10*max(abs(eig(b)))*max(diag(abs(d)))*feps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a,b);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t1all = t1all & t1 & tb;
//end
//% Now the double version
//t2all = 1;
//for i=2:4:100
// a = double(randn(i)); a = a + a';
// b = double(randn(i)); b = b*b';
// [v,d] = eig(a,b);
// emat = a*v - b*v*d;
// er = max(abs(emat(:)));
// bnd = 10*max(abs(eig(b)))*max(diag(abs(d)))*eps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a,b);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t2all = t2all & t1 & tb;
//end
//% Now the complex version
//t3all = 1;
//for i=2:4:100
// a = complex(randn(i)+j*randn(i)); a = a + a';
// b = complex(randn(i)+j*randn(i)); b = b*b';
// [v,d] = eig(a,b);
// emat = a*v - b*v*d;
// er = max(abs(emat(:)));
// bnd = 10*max(abs(eig(b)))*max(diag(abs(d)))*feps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a,b);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t3all = t3all & t1 & tb;
//end
//% Now the double version
//t4all = 1;
//for i=2:4:100
// a = dcomplex(randn(i)+j*randn(i)); a = a + a';
// b = dcomplex(randn(i)+j*randn(i)); b = b*b';
// [v,d] = eig(a,b);
// emat = a*v - b*v*d;
// er = max(abs(emat(:)));
// bnd = 10*max(abs(eig(b)))*max(diag(abs(d)))*eps*i;
// t1 = (er < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',er,bnd,i); end
// g = eig(a,b);
// e2 = max(abs(sort(g)-sort(diag(d))));
// tb = e2<bnd;
// if (~tb) printf('compact/full decomp mismatch: er = %e (num = %d)\n',e2,i); end
// t4all = t4all & t1 & tb;
//end
//t = t1all & t2all & t3all & t4all;
//@}
//@{ test_eig6.m
//function t = test_eig6
// try
// eig([1,2;3,4],eye(1));
// catch
// end
// t = true;
//@}
//!
ArrayVector EigFunction(int nargout, const ArrayVector& arg, Interpreter* m_eval) {
bool balance;
if (arg.size() == 0)
throw Exception("eig function requires at least one argument");
if (arg.size() == 1)
balance = true;
else {
Array b(arg[1]);
if (b.isString()) {
string b2 = b.getContentsAsStringUpper();
if (b2=="NOBALANCE")
balance = false;
}
else
return GenEigFunction(nargout, arg, m_eval);
}
Array A(arg[0]);
if (!A.is2D())
throw Exception("Cannot apply matrix operations to N-Dimensional arrays.");
if (A.anyNotFinite())
throw Exception("eig only defined for matrices with finite entries.");
ArrayVector retval;
Array V, D;
if (nargout > 1) {
if (A.isSymmetric())
EigenDecomposeFullSymmetric(A,V,D, m_eval);
else
EigenDecomposeFullGeneral(A,V,D,balance, m_eval);
retval.push_back(V);
retval.push_back(D);
} else {
if (A.isSymmetric())
EigenDecomposeCompactSymmetric(A,D, m_eval);
else
EigenDecomposeCompactGeneral(A,D,balance, m_eval);
retval.push_back(D);
}
return retval;
}
//!
//@Module EIGS Sparse Matrix Eigendecomposition
//@@Section SPARSE
//@@Usage
//Computes the eigendecomsition of a sparse square matrix. The
//@|eigs| function has several forms. The most general form is
//@[
// [V,D] = eigs(A,k,sigma)
//@]
//where @|A| is the matrix to analyze, @|k| is the number of
//eigenvalues to compute and @|sigma| determines which eigenvallues
//to solve for. Valid values for @|sigma| are
// 'lm' - largest magnitude
// 'sm' - smallest magnitude
// 'la' - largest algebraic (for real symmetric problems)
// 'sa' - smallest algebraic (for real symmetric problems)
// 'be' - both ends (for real symmetric problems)
// 'lr' - largest real part
// 'sr' - smallest real part
// 'li' - largest imaginary part
// 'si' - smallest imaginary part
// scalar - find the eigenvalues closest to @|sigma|.
//The returned matrix @|V| contains the eigenvectors, and @|D|
//stores the eigenvalues. The related form
//@[
// d = eigs(A,k,sigma)
//@]
//computes only the eigenvalues and not the eigenvectors. If @|sigma|
//is omitted, as in the forms
//@[
// [V,D] = eigs(A,k)
//@]
//and
//@[
// d = eigs(A,k)
//@]
//then @|eigs| returns the largest magnitude eigenvalues (and optionally
//the associated eigenvectors). As an even simpler form, the forms
//@[
// [V,D] = eigs(A)
//@]
//and
//@[
// d = eigs(A)
//@]
//then @|eigs| returns the six largest magnitude eigenvalues of @|A| and
//optionally the eigenvectors. The @|eigs| function uses ARPACK to
//compute the eigenvectors and/or eigenvalues. Note that due to a
//limitation in the interface into ARPACK from FreeMat, the number of
//eigenvalues that are to be computed must be strictly smaller than the
//number of columns (or rows) in the matrix.
//@@Example
//Here is an example of using @|eigs| to calculate eigenvalues
//of a matrix, and a comparison of the results with @|eig|
//@<
//a = sparse(rand(9))
//eigs(a)
//eig(full(a))
//@>
//Next, we exercise some of the variants of @|eigs|:
//@<
//eigs(a,4,'sm')
//eigs(a,4,'lr')
//eigs(a,4,'sr')
//@>
//@@Tests
//@{ test_sparse45.m
//function x = test_sparse45
//a = rand(10);
//A = sparse(a);
//p = eigs(A,4,0.634);
//d = eig(a);
//e = min(abs(d*ones(1,4) - ones(10,1)*p.'));
//x = max(e) < 100*eps;
//@}
//!
ArrayVector EigsFunction(int nargout, const ArrayVector& arg) {
if (arg.size() == 0)
throw Exception("eigs function requires at least one argument");
Array A(arg[0]);
if (!A.sparse())
throw Exception("eigs only applies to sparse matrix arguments");
int k;
if (A.getDimensionLength(0) != A.getDimensionLength(1))
throw Exception("eigs can only be applied to square matrices.");
if (arg.size() < 2) {
k = 6;
if (k >= A.getDimensionLength(0))
k = A.getDimensionLength(0) - 1;
} else {
Array kval(arg[1]);
k = kval.getContentsAsIntegerScalar();
}
if (A.isComplex())
A.promoteType(FM_DCOMPLEX);
else
A.promoteType(FM_DOUBLE);
bool shiftFlag;
string whichflag;
double sigma[2];
if (arg.size() < 3) {
shiftFlag = false;
whichflag = "LM";
} else {
Array S(arg[2]);
if (S.isString()) {
shiftFlag = false;
string stxt = S.getContentsAsStringUpper();
if ((stxt == "LM") || (stxt == "SM") || (stxt == "LA") || (stxt == "SA") ||
(stxt == "BE") || (stxt == "LR") || (stxt == "SR") || (stxt == "LI") ||
(stxt == "SI"))
whichflag = stxt;
else
throw Exception("Unrecognized option for sigma - it must be either 'lm', 'sm', 'la', 'sa', 'be', 'lr', 'sr', 'li', or 'si'");
} else {
if (!S.isScalar())
throw Exception("shift parameter sigma must be a scalar");
if (S.isComplex()) {
S.promoteType(FM_DCOMPLEX);
double *cp = (double*) S.getDataPointer();
sigma[0] = cp[0];
sigma[1] = cp[1];
} else {
sigma[0] = S.getContentsAsDoubleScalar();
sigma[1] = 0;
}
shiftFlag = true;
}
}
if (!shiftFlag)
return SparseEigDecompose(nargout,A,k,whichflag);
else
return SparseEigDecomposeShifted(nargout,A,k,sigma);
}
ArrayVector QRDNoPivotFunction(bool compactDec, Array A) {
int nrows = A.getDimensionLength(0);
int ncols = A.getDimensionLength(1);
int orows = nrows;
int ocols = ncols;
if ((!compactDec) && (nrows > ncols)) {
Dimensions newDim(nrows,nrows);
A.resize(newDim);
ncols = nrows;
} else
compactDec = true;
Class Aclass(A.dataClass());
if (Aclass < FM_FLOAT) {
A.promoteType(FM_DOUBLE);
Aclass = FM_DOUBLE;
}
int minmn = (nrows < ncols) ? nrows : ncols;
ArrayVector retval;
Array rmat, qmat;
Dimensions dim;
switch (Aclass) {
case FM_FLOAT:
{
float *q = (float*) Malloc(nrows*minmn*sizeof(float));
float *r = (float*) Malloc(ncols*minmn*sizeof(float));
floatQRD(nrows,ncols,q,r,(float*) A.getReadWriteDataPointer());
if (!compactDec) {
float *r2 = (float*) Malloc(orows*ocols*sizeof(float));
memcpy(r2,r,orows*ocols*sizeof(float));
dim.set(0,orows);
dim.set(1,ocols);
rmat = Array(FM_FLOAT,dim,r2);
Free(r);
} else {
dim.set(0,minmn);
dim.set(1,ncols);
rmat = Array(FM_FLOAT,dim,r);
}
dim.set(0,nrows);
dim.set(1,minmn);
qmat = Array(FM_FLOAT,dim,q);
retval.push_back(qmat);
retval.push_back(rmat);
break;
}
case FM_DOUBLE:
{
double *q = (double*) Malloc(nrows*minmn*sizeof(double));
double *r = (double*) Malloc(ncols*minmn*sizeof(double));
doubleQRD(nrows,ncols,q,r,(double*) A.getReadWriteDataPointer());
if (!compactDec) {
double *r2 = (double*) Malloc(orows*ocols*sizeof(double));
memcpy(r2,r,orows*ocols*sizeof(double));
dim.set(0,orows);
dim.set(1,ocols);
rmat = Array(FM_DOUBLE,dim,r2);
Free(r);
} else {
dim.set(0,minmn);
dim.set(1,ncols);
rmat = Array(FM_DOUBLE,dim,r);
}
dim.set(0,nrows);
dim.set(1,minmn);
qmat = Array(FM_DOUBLE,dim,q);
retval.push_back(qmat);
retval.push_back(rmat);
break;
}
case FM_COMPLEX:
{
float *q = (float*) Malloc(2*nrows*minmn*sizeof(float));
float *r = (float*) Malloc(2*ncols*minmn*sizeof(float));
complexQRD(nrows,ncols,q,r,(float*) A.getReadWriteDataPointer());
if (!compactDec) {
float *r2 = (float*) Malloc(2*orows*ocols*sizeof(float));
memcpy(r2,r,2*orows*ocols*sizeof(float));
dim.set(0,orows);
dim.set(1,ocols);
rmat = Array(FM_COMPLEX,dim,r2);
Free(r);
} else {
dim.set(0,minmn);
dim.set(1,ncols);
rmat = Array(FM_COMPLEX,dim,r);
}
dim.set(0,nrows);
dim.set(1,minmn);
qmat = Array(FM_COMPLEX,dim,q);
retval.push_back(qmat);
retval.push_back(rmat);
break;
}
case FM_DCOMPLEX:
{
double *q = (double*) Malloc(2*nrows*minmn*sizeof(double));
double *r = (double*) Malloc(2*ncols*minmn*sizeof(double));
dcomplexQRD(nrows,ncols,q,r,(double*) A.getReadWriteDataPointer());
if (!compactDec) {
double *r2 = (double*) Malloc(2*orows*ocols*sizeof(double));
memcpy(r2,r,2*orows*ocols*sizeof(double));
dim.set(0,orows);
dim.set(1,ocols);
rmat = Array(FM_DCOMPLEX,dim,r2);
Free(r);
} else {
dim.set(0,minmn);
dim.set(1,ncols);
rmat = Array(FM_DCOMPLEX,dim,r);
}
dim.set(0,nrows);
dim.set(1,minmn);
qmat = Array(FM_DCOMPLEX,dim,q);
retval.push_back(qmat);
retval.push_back(rmat);
break;
}
}
return retval;
}
ArrayVector QRDPivotFunction(bool compactDec, Array A) {
int nrows = A.getDimensionLength(0);
int ncols = A.getDimensionLength(1);
int orows = nrows;
int ocols = ncols;
int i;
bool compactSav = compactDec;
if ((!compactDec) && (nrows > ncols)) {
Dimensions newDim(nrows,nrows);
A.resize(newDim);
ncols = nrows;
} else
compactDec = true;
Class Aclass(A.dataClass());
if (Aclass < FM_FLOAT) {
A.promoteType(FM_DOUBLE);
Aclass = FM_DOUBLE;
}
int minmn = (nrows < ncols) ? nrows : ncols;
ArrayVector retval;
Array rmat, qmat, pmat;
Dimensions dim;
switch (Aclass) {
case FM_FLOAT:
{
float *q = (float*) Malloc(nrows*minmn*sizeof(float));
float *r = (float*) Malloc(ncols*minmn*sizeof(float));
int *p = (int*) Malloc(ncols*sizeof(int));
floatQRDP(nrows,ncols,q,r,p,(float*) A.getReadWriteDataPointer());
if (!compactDec) {
float *r2 = (float*) Malloc(orows*ocols*sizeof(float));
memcpy(r2,r,orows*ocols*sizeof(float));
dim.set(0,orows);
dim.set(1,ocols);
rmat = Array(FM_FLOAT,dim,r2);
Free(r);
} else {
dim.set(0,minmn);
dim.set(1,ncols);
rmat = Array(FM_FLOAT,dim,r);
}
if (!compactSav) {
int *p2 = (int*) Malloc(ncols*ncols*sizeof(int));
for (i=0;i<ncols;i++)
p2[p[i] + i*ncols - 1] = 1;
dim.set(0,ncols);
dim.set(1,ncols);
pmat = Array(FM_INT32,dim,p2);
Free(p);
} else {
dim.set(0,1);
dim.set(1,ncols);
pmat = Array(FM_INT32,dim,p);
}
dim.set(0,nrows);
dim.set(1,minmn);
qmat = Array(FM_FLOAT,dim,q);
retval.push_back(qmat);
retval.push_back(rmat);
retval.push_back(pmat);
break;
}
case FM_DOUBLE:
{
double *q = (double*) Malloc(nrows*minmn*sizeof(double));
double *r = (double*) Malloc(ncols*minmn*sizeof(double));
int *p = (int*) Malloc(ncols*sizeof(int));
doubleQRDP(nrows,ncols,q,r,p,(double*) A.getReadWriteDataPointer());
if (!compactDec) {
double *r2 = (double*) Malloc(orows*ocols*sizeof(double));
memcpy(r2,r,orows*ocols*sizeof(double));
dim.set(0,orows);
dim.set(1,ocols);
rmat = Array(FM_DOUBLE,dim,r2);
Free(r);
} else {
dim.set(0,minmn);
dim.set(1,ncols);
rmat = Array(FM_DOUBLE,dim,r);
}
if (!compactSav) {
int *p2 = (int*) Malloc(ncols*ncols*sizeof(int));
for (i=0;i<ncols;i++)
p2[p[i] + i*ncols - 1] = 1;
dim.set(0,ncols);
dim.set(1,ncols);
pmat = Array(FM_INT32,dim,p2);
Free(p);
} else {
dim.set(0,1);
dim.set(1,ncols);
pmat = Array(FM_INT32,dim,p);
}
dim.set(0,nrows);
dim.set(1,minmn);
qmat = Array(FM_DOUBLE,dim,q);
retval.push_back(qmat);
retval.push_back(rmat);
retval.push_back(pmat);
break;
}
case FM_COMPLEX:
{
float *q = (float*) Malloc(2*nrows*minmn*sizeof(float));
float *r = (float*) Malloc(2*ncols*minmn*sizeof(float));
int *p = (int*) Malloc(ncols*sizeof(int));
complexQRDP(nrows,ncols,q,r,p,(float*) A.getReadWriteDataPointer());
if (!compactDec) {
float *r2 = (float*) Malloc(2*orows*ocols*sizeof(float));
memcpy(r2,r,2*orows*ocols*sizeof(float));
dim.set(0,orows);
dim.set(1,ocols);
rmat = Array(FM_COMPLEX,dim,r2);
Free(r);
} else {
dim.set(0,minmn);
dim.set(1,ncols);
rmat = Array(FM_COMPLEX,dim,r);
}
if (!compactSav) {
int *p2 = (int*) Malloc(ncols*ncols*sizeof(int));
for (i=0;i<ncols;i++)
p2[p[i] + i*ncols - 1] = 1;
dim.set(0,ncols);
dim.set(1,ncols);
pmat = Array(FM_INT32,dim,p2);
Free(p);
} else {
dim.set(0,1);
dim.set(1,ncols);
pmat = Array(FM_INT32,dim,p);
}
dim.set(0,nrows);
dim.set(1,minmn);
qmat = Array(FM_COMPLEX,dim,q);
retval.push_back(qmat);
retval.push_back(rmat);
retval.push_back(pmat);
break;
}
case FM_DCOMPLEX:
{
double *q = (double*) Malloc(2*nrows*minmn*sizeof(double));
double *r = (double*) Malloc(2*ncols*minmn*sizeof(double));
int *p = (int*) Malloc(ncols*sizeof(int));
dcomplexQRDP(nrows,ncols,q,r,p,(double*) A.getReadWriteDataPointer());
if (!compactDec) {
double *r2 = (double*) Malloc(2*orows*ocols*sizeof(double));
memcpy(r2,r,2*orows*ocols*sizeof(double));
dim.set(0,orows);
dim.set(1,ocols);
rmat = Array(FM_DCOMPLEX,dim,r2);
Free(r);
} else {
dim.set(0,minmn);
dim.set(1,ncols);
rmat = Array(FM_DCOMPLEX,dim,r);
}
if (!compactSav) {
int *p2 = (int*) Malloc(ncols*ncols*sizeof(int));
for (i=0;i<ncols;i++)
p2[p[i] + i*ncols - 1] = 1;
dim.set(0,ncols);
dim.set(1,ncols);
pmat = Array(FM_INT32,dim,p2);
Free(p);
} else {
dim.set(0,1);
dim.set(1,ncols);
pmat = Array(FM_INT32,dim,p);
}
dim.set(0,nrows);
dim.set(1,minmn);
qmat = Array(FM_DCOMPLEX,dim,q);
retval.push_back(qmat);
retval.push_back(rmat);
retval.push_back(pmat);
break;
}
}
return retval;
}
//!
//@Module QR QR Decomposition of a Matrix
//@@Section TRANSFORMS
//@@Usage
//Computes the QR factorization of a matrix. The @|qr| function has
//multiple forms, with and without pivoting. The non-pivot version
//has two forms, a compact version and a full-blown decomposition
//version. The compact version of the decomposition of a matrix
//of size @|M x N| is
//@[
// [q,r] = qr(a,0)
//@]
//where @|q| is a matrix of size @|M x L| and @|r| is a matrix of
//size @|L x N| and @|L = min(N,M)|, and @|q*r = a|. The QR decomposition is
//such that the columns of @|Q| are orthonormal, and @|R| is upper
//triangular. The decomposition is computed using the LAPACK
//routine @|xgeqrf|, where @|x| is the precision of the matrix. Unlike
//MATLAB (and other MATLAB-compatibles), FreeMat supports decompositions
//of all four floating point types, @|float, complex, double, dcomplex|.
//
//The second form of the non-pivot decomposition omits the second @|0|
//argument:
//@[
// [q,r] = qr(a)
//@]
//This second form differs from the previous form only for matrices
//with more rows than columns (@|M > N|). For these matrices, the
//full decomposition is of a matrix @|Q| of size @|M x M| and
//a matrix @|R| of size @|M x N|. The full decomposition is computed
//using the same LAPACK routines as the compact decomposition, but
//on an augmented matrix @|[a 0]|, where enough columns are added to
//form a square matrix.
//
//Generally, the QR decomposition will not return a matrix @|R| with
//diagonal elements in any specific order. The remaining two forms
//of the @|qr| command utilize permutations of the columns of @|a|
//so that the diagonal elements of @|r| are in decreasing magnitude.
//To trigger this form of the decomposition, a third argument is
//required, which records the permutation applied to the argument @|a|.
//The compact version is
//@[
// [q,r,e] = qr(a,0)
//@]
//where @|e| is an integer vector that describes the permutation of
//the columns of @|a| necessary to reorder the diagonal elements of
//@|r|. This result is computed using the LAPACK routines @|(s,d)geqp3|.
//In the non-compact version of the QR decomposition with pivoting,
//@[
// [q,r,e] = qr(a)
//@]
//the returned matrix @|e| is a permutation matrix, such that
//@|q*r*e' = a|.
//@@Tests
//@{ test_qr1.m
//% Test the QR decomposition (without pivoting, minimal decomposition, double)
//function test_val = test_qr1
//a = [1,2,3;4,5,6;7,8,9;10,0,5];
//% double precision
//[q,r] = qr(a,0);
//qorg = q'*q - diag(ones(1,3));
//t1 = (max(abs(qorg(:)))<1e-15);
//t2 = ((size(q,1) == 4) & (size(q,2) == 3));
//t3 = ((size(r,1) == 3) & (size(r,2) == 3));
//aorg = q*r - a;
//t4 = (max(abs(qorg(:)))<1e-15);
//t5 = max(diag(r,-1)) == 0;
//t6 = max(diag(r,-2)) == 0;
//test_val = test(t1) & test(t2) & test(t3) & test(t4) & test(t5) & test(t6);
//@}
//@{ test_qr2.m
//% Test the QR decomposition (without pivoting, minimal decomposition, float)
//function test_val = test_qr2
//a = float([1,2,3;4,5,6;7,8,9;10,0,5]);
//% float precision
//[q,r] = qr(a,0);
//qorg = q'*q - diag(ones(1,3));
//t1 = (max(abs(qorg(:)))<1e-6);
//t2 = ((size(q,1) == 4) & (size(q,2) == 3));
//t3 = ((size(r,1) == 3) & (size(r,2) == 3));
//aorg = q*r - a;
//t4 = (max(abs(qorg(:)))<1e-6);
//t5 = max(diag(r,-1)) == 0;
//t6 = max(diag(r,-2)) == 0;
//test_val = test(t1) & test(t2) & test(t3) & test(t4) & test(t5) & test(t6);
//@}
//@{ test_qr3.m
//% Test the QR decomposition (without pivoting, minimal decomposition, complex)
//function test_val = test_qr3
//a = float([1,2,3;4,5,6;7,8,9;10,0,5]);
//a = a + i*(a-1);
//% float precision
//[q,r] = qr(a,0);
//qorg = q'*q - diag(ones(1,3));
//t1 = (max(abs(qorg(:)))<1e-6);
//t2 = ((size(q,1) == 4) & (size(q,2) == 3));
//t3 = ((size(r,1) == 3) & (size(r,2) == 3));
//aorg = q*r - a;
//t4 = (max(abs(qorg(:)))<1e-6);
//t5 = max(diag(r,-1)) == 0;
//t6 = max(diag(r,-2)) == 0;
//test_val = test(t1) & test(t2) & test(t3) & test(t4) & test(t5) & test(t6);
//@}
//@{ test_qr4.m
//% Test the QR decomposition (without pivoting, minimal decomposition, dcomplex)
//function test_val = test_qr4
//a = double([1,2,3;4,5,6;7,8,9;10,0,5]);
//a = a + i*(a-1);
//% float precision
//[q,r] = qr(a,0);
//qorg = q'*q - diag(ones(1,3));
//t1 = (max(abs(qorg(:)))<1e-15);
//t2 = ((size(q,1) == 4) & (size(q,2) == 3));
//t3 = ((size(r,1) == 3) & (size(r,2) == 3));
//aorg = q*r - a;
//t4 = (max(abs(qorg(:)))<1e-15);
//t5 = max(diag(r,-1)) == 0;
//t6 = max(diag(r,-2)) == 0;
//test_val = test(t1) & test(t2) & test(t3) & test(t4) & test(t5) & test(t6);
//@}
//@{ test_qr5.m
//% Test the QR decomposition (without pivoting, double)
//function test_val = test_qr5
//a = [1,2,3;4,5,6;7,8,9;10,0,5];
//% double precision
//[q,r] = qr(a);
//qorg = q'*q - diag(ones(1,4));
//t1 = (max(abs(qorg(:)))<1e-15);
//t2 = ((size(q,1) == 4) & (size(q,2) == 4));
//t3 = ((size(r,1) == 4) & (size(r,2) == 3));
//aorg = q*r - a;
//t4 = (max(abs(qorg(:)))<1e-15);
//t5 = max(diag(r,-1)) == 0;
//t6 = max(diag(r,-2)) == 0;
//t7 = max(diag(r,-3)) == 0;
//test_val = test(t1) & test(t2) & test(t3) & test(t4) & test(t5) & test(t6) & test(t7);
//@}
//@{ test_qr6.m
//% Test the QR decomposition (without pivoting, float)
//function test_val = test_qr6
//a = float([1,2,3;4,5,6;7,8,9;10,0,5]);
//% float precision
//[q,r] = qr(a);
//qorg = q'*q - diag(ones(1,4));
//t1 = (max(abs(qorg(:)))<1e-6);
//t2 = ((size(q,1) == 4) & (size(q,2) == 4));
//t3 = ((size(r,1) == 4) & (size(r,2) == 3));
//aorg = q*r - a;
//t4 = (max(abs(qorg(:)))<1e-6);
//t5 = max(diag(r,-1)) == 0;
//t6 = max(diag(r,-2)) == 0;
//t7 = max(diag(r,-3)) == 0;
//test_val = test(t1) & test(t2) & test(t3) & test(t4) & test(t5) & test(t6) & test(t7);
//@}
//@{ test_qr7.m
//% Test the QR decomposition (without pivoting, complex)
//function test_val = test_qr7
//a = float([1,2,3;4,5,6;7,8,9;10,0,5]);
//a = a + i*(a-1);
//% float precision
//[q,r] = qr(a);
//qorg = q'*q - diag(ones(1,4));
//t1 = (max(abs(qorg(:)))<1e-6);
//t2 = ((size(q,1) == 4) & (size(q,2) == 4));
//t3 = ((size(r,1) == 4) & (size(r,2) == 3));
//aorg = q*r - a;
//t4 = (max(abs(qorg(:)))<1e-6);
//t5 = max(diag(r,-1)) == 0;
//t6 = max(diag(r,-2)) == 0;
//t7 = max(diag(r,-3)) == 0;
//test_val = test(t1) & test(t2) & test(t3) & test(t4) & test(t5) & test(t6) & test(t7);
//@}
//@{ test_qr8.m
//% Test the QR decomposition (without pivoting, dcomplex)
//function test_val = test_qr8
//a = double([1,2,3;4,5,6;7,8,9;10,0,5]);
//a = a + i*(a-1);
//% float precision
//[q,r] = qr(a);
//qorg = q'*q - diag(ones(1,4));
//t1 = (max(abs(qorg(:)))<1e-15);
//t2 = ((size(q,1) == 4) & (size(q,2) == 4));
//t3 = ((size(r,1) == 4) & (size(r,2) == 3));
//aorg = q*r - a;
//t4 = (max(abs(qorg(:)))<1e-15);
//t5 = max(diag(r,-1)) == 0;
//t6 = max(diag(r,-2)) == 0;
//t7 = max(diag(r,-3)) == 0;
//test_val = test(t1) & test(t2) & test(t3) & test(t4) & test(t5) & test(t6) & test(t7);
//@}
//@{ test_qr9.m
//% Test the QR decomposition (full, without pivoting)
//function test_val = test_qr9
// t1all = 1;
// for n=2:4:100
// a = float(rand(n,2*n));
// [q,r] = qr(a);
// err = abs(a-q*r);
// err = max(err(:));
// bnd = 8*max(abs(diag(r)))*feps*(n^(0.5));
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// t2all = 1;
// for n=2:4:100
// a = double(rand(n,2*n));
// [q,r] = qr(a);
// err = abs(a-q*r);
// err = max(err(:));
// bnd = 8*max(abs(diag(r)))*eps*(n^(0.5));
// t2 = (err < bnd);
// if (~t2) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t2all = t2all & t2;
// end
// t3all = 1;
// for n=2:4:100
// a = complex(rand(n,2*n)+i*rand(n,2*n));
// [q,r] = qr(a);
// err = abs(a-q*r);
// err = max(err(:));
// bnd = 8*max(abs(diag(r)))*feps*(n^(0.5));
// t3 = (err < bnd);
// if (~t3) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t3all = t3all & t3;
// end
// t4all = 1;
// for n=2:4:100
// a = dcomplex(rand(n,2*n)+i*rand(n,2*n));
// [q,r] = qr(a);
// err = abs(a-q*r);
// err = max(err(:));
// bnd = 8*max(abs(diag(r)))*eps*(n^(0.5));
// t4 = (err < bnd);
// if (~t4) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t4all = t4all & t4;
// end
// test_val = t1all & t2all & t3all & t4all;
//@}
//@{ test_qr10.m
//% Test the QR decomposition (compact, without pivoting)
//function test_val = test_qr10
// t1all = 1;
// for n=2:4:100
// a = float(rand(2*n,n));
// [q,r] = qr(a,0);
// err = abs(a-q*r);
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*feps*(n^(0.5));
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// t2all = 1;
// for n=2:4:100
// a = double(rand(2*n,n));
// [q,r] = qr(a,0);
// err = abs(a-q*r);
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*eps*(n^(0.5));
// t2 = (err < bnd);
// if (~t2) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t2all = t2all & t2;
// end
// t3all = 1;
// for n=2:4:100
// a = complex(rand(2*n,n)+i*rand(2*n,n));
// [q,r] = qr(a,0);
// err = abs(a-q*r);
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*feps*(n^(0.5));
// t3 = (err < bnd);
// if (~t3) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t3all = t3all & t3;
// end
// t4all = 1;
// for n=2:4:100
// a = dcomplex(rand(2*n,n)+i*rand(2*n,n));
// [q,r] = qr(a,0);
// err = abs(a-q*r);
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*eps*(n^(0.5));
// t4 = (err < bnd);
// if (~t4) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t4all = t4all & t4;
// end
// test_val = t1all & t2all & t3all & t4all;
//@}
//@{ test_qr11.m
//% Test the QR decomposition (compact, pivoting)
//function test_val = test_qr11
// t1all = 1;
// for n=2:4:100
// a = float(rand(2*n,n));
// [q,r,ec] = qr(a,0);
// e = zeros(n,n);
// p = 0:(n-1);
// mdx = ec + p*n;
// e(mdx) = 1;
// err = abs(a-q*r*e');
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*feps*(n^(0.5));
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// t2all = 1;
// for n=2:4:100
// a = double(rand(2*n,n));
// [q,r,ec] = qr(a,0);
// e = zeros(n,n);
// p = 0:(n-1);
// mdx = ec + p*n;
// e(mdx) = 1;
// err = abs(a-q*r*e');
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*eps*(n^(0.5));
// t2 = (err < bnd);
// if (~t2) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t2all = t2all & t2;
// end
// t3all = 1;
// for n=2:4:100
// a = complex(rand(2*n,n)+i*rand(2*n,n));
// [q,r,ec] = qr(a,0);
// e = zeros(n,n);
// p = 0:(n-1);
// mdx = ec + p*n;
// e(mdx) = 1;
// err = abs(a-q*r*e');
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*feps*(n^(0.5));
// t3 = (err < bnd);
// if (~t3) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t3all = t3all & t3;
// end
// t4all = 1;
// for n=2:4:100
// a = dcomplex(rand(2*n,n)+i*rand(2*n,n));
// [q,r,ec] = qr(a,0);
// e = zeros(n,n);
// p = 0:(n-1);
// mdx = ec + p*n;
// e(mdx) = 1;
// err = abs(a-q*r*e');
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*eps*(n^(0.5));
// t4 = (err < bnd);
// if (~t4) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t4all = t4all & t4;
// end
// test_val = t1all & t2all & t3all & t4all;
//@}
//@{ test_qr12.m
//% Test the QR decomposition (full, pivoting)
//function test_val = test_qr12
// t1all = 1;
// for n=2:4:100
// a = float(rand(n,n));
// [q,r,e] = qr(a);
// err = abs(a-q*r*e');
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*feps*(n^(0.5));
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// t2all = 1;
// for n=2:4:100
// a = double(rand(n,n));
// [q,r,e] = qr(a);
// err = abs(a-q*r*e');
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*eps*(n^(0.5));
// t2 = (err < bnd);
// if (~t2) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t2all = t2all & t2;
// end
// t3all = 1;
// for n=2:4:100
// a = complex(rand(n,n)+i*rand(n,n));
// [q,r,e] = qr(a);
// err = abs(a-q*r*e');
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*feps*(n^(0.5));
// t3 = (err < bnd);
// if (~t3) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t3all = t3all & t3;
// end
// t4all = 1;
// for n=2:4:100
// a = dcomplex(rand(n,n)+i*rand(n,n));
// [q,r,e] = qr(a);
// err = abs(a-q*r*e');
// err = max(err(:));
// bnd = 10*max(abs(diag(r)))*eps*(n^(0.5));
// t4 = (err < bnd);
// if (~t4) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t4all = t4all & t4;
// end
// test_val = t1all & t2all & t3all & t4all;
//@}
//!
ArrayVector QRDFunction(int nargout, const ArrayVector& arg) {
if (arg.size() < 1)
throw Exception("qr function requires at least one argument - the matrix to decompose.");
Array A(arg[0]);
// Test for numeric
if (A.isReferenceType())
throw Exception("Cannot apply qr decomposition to reference types.");
if (!A.is2D())
throw Exception("Cannot apply matrix operations to N-Dimensional arrays.");
if (A.anyNotFinite())
throw Exception("QR Decomposition only defined for matrices with finite entries.");
int nrows = A.getDimensionLength(0);
int ncols = A.getDimensionLength(1);
bool compactDecomposition = false;
if (arg.size() == 2) {
Array cflag(arg[1]);
int cflag_int = cflag.getContentsAsIntegerScalar();
if (cflag_int == 0) compactDecomposition = true;
}
if (nargout == 3)
return QRDPivotFunction(compactDecomposition, A);
else
return QRDNoPivotFunction(compactDecomposition, A);
}
//!
//@Module SVD Singular Value Decomposition of a Matrix
//@@Section TRANSFORMS
//@@Usage
//Computes the singular value decomposition (SVD) of a matrix. The
//@|svd| function has three forms. The first returns only the singular
//values of the matrix:
//@[
// s = svd(A)
//@]
//The second form returns both the singular values in a diagonal
//matrix @|S|, as well as the left and right eigenvectors.
//@[
// [U,S,V] = svd(A)
//@]
//The third form returns a more compact decomposition, with the
//left and right singular vectors corresponding to zero singular
//values being eliminated. The syntax is
//@[
// [U,S,V] = svd(A,0)
//@]
//@@Function Internals
//Recall that @|sigma_i| is a singular value of an @|M x N|
//matrix @|A| if there exists two vectors @|u_i, v_i| where @|u_i| is
//of length @|M|, and @|v_i| is of length @|u_i| and
//\[
// A v_i = \sigma_i u_i
//\]
//and generally
//\[
// A = \sum_{i=1}^{K} \sigma_i u_i*v_i',
//\]
//where @|K| is the rank of @|A|. In matrix form, the left singular
//vectors @|u_i| are stored in the matrix @|U| as
//\[
// U = [u_1,\ldots,u_m], V = [v_1,\ldots,v_n]
//\]
//The matrix @|S| is then of size @|M x N| with the singular
//values along the diagonal. The SVD is computed using the
//@|LAPACK| class of functions @|GESDD|.
//@@Examples
//Here is an example of a partial and complete singular value
//decomposition.
//@<
//A = float(randn(2,3))
//[U,S,V] = svd(A)
//U*S*V'
//svd(A)
//@>
//@@Tests
//@{ test_svd1.m
//% Test the svd with finite and infinite args (functional only)
//function test_val = test_svd1
// a = [1,3;4,2];
// [u,s,v] = svd(a);
// res1 = 1;
// res2 = 0;
// res3 = 0;
// res4 = 0;
// try
// c = a;
// c(2,2) = inf;
// svd(c);
// catch
// res2 = 1;
// end
// try
// c = a;
// c(2,2) = -inf;
// svd(c);
// catch
// res3 = 1;
// end
// try
// c = a;
// c(2,2) = nan;
// svd(c);
// catch
// res4 = 1;
// end
//
// test_val = res1 & res2 & res3 & res4;
//@}
//@{ test_svd2.m
//% Test the svd with finite and infinite args (functional only)
//function test_val = test_svd2
// a = double([1,3;4,2]);
// [u,s,v] = svd(a);
// res1 = 1;
// res2 = 0;
// res3 = 0;
// res4 = 0;
// try
// c = a;
// c(2,2) = inf;
// svd(c);
// catch
// res2 = 1;
// end
// try
// c = a;
// c(2,2) = -inf;
// svd(c);
// catch
// res3 = 1;
// end
// try
// c = a;
// c(2,2) = nan;
// svd(c);
// catch
// res4 = 1;
// end
//
// test_val = res1 & res2 & res3 & res4;
//@}
//@{ test_svd3.m
//% Test the svd with finite and infinite args (functional only)
//function test_val = test_svd3
// a = complex([1,3;4,2]);
// [u,s,v] = svd(a);
// res1 = 1;
// res2 = 0;
// res3 = 0;
// res4 = 0;
// try
// c = a;
// c(2,2) = inf;
// svd(c);
// catch
// res2 = 1;
// end
// try
// c = a;
// c(2,2) = -inf;
// svd(c);
// catch
// res3 = 1;
// end
// try
// c = a;
// c(2,2) = nan;
// svd(c);
// catch
// res4 = 1;
// end
//
// test_val = res1 & res2 & res3 & res4;
//@}
//@{ test_svd4.m
//% Test the svd with finite and infinite args (functional only)
//function test_val = test_svd4
// a = dcomplex([1,3;4,2]);
// [u,s,v] = svd(a);
// res1 = 1;
// res2 = 0;
// res3 = 0;
// res4 = 0;
// try
// c = a;
// c(2,2) = inf;
// svd(c);
// catch
// res2 = 1;
// end
// try
// c = a;
// c(2,2) = -inf;
// svd(c);
// catch
// res3 = 1;
// end
// try
// c = a;
// c(2,2) = nan;
// svd(c);
// catch
// res4 = 1;
// end
//
// test_val = res1 & res2 & res3 & res4;
//@}
//@{ test_svd5.m
//% Test the full version of the svd with float vectors
//function test_val = test_svd5
// t1all = 1;
// p = [];
// for n=2:4:100
// a = float(rand(n,2*n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// p(n) = err;
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = float(rand(2*n,n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = float(rand(n,n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// test_val = t1all;
//@}
//@{ test_svd6.m
//% Test the full version of the svd with complex vectors
//function test_val = test_svd6
// t1all = 1;
// p = [];
// for n=2:4:100
// a = float(rand(n,2*n))+i*float(rand(n,2*n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// p(n) = err;
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = float(rand(2*n,n))+i*float(rand(2*n,n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = float(rand(n,n))+i*float(rand(n,n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// test_val = t1all;
//@}
//@{ test_svd7.m
//% Test the full version of the svd with double vectors
//function test_val = test_svd7
// t1all = 1;
// p = [];
// for n=2:4:100
// a = double(rand(n,2*n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// p(n) = err;
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = double(rand(2*n,n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = double(rand(n,n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// test_val = t1all;
//@}
//@{ test_svd8.m
//% Test the full version of the svd with dcomplex vectors
//function test_val = test_svd8
// t1all = 1;
// p = [];
// for n=2:4:100
// a = double(rand(n,2*n))+i*double(rand(n,2*n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// p(n) = err;
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = double(rand(2*n,n))+i*double(rand(2*n,n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = double(rand(n,n))+i*double(rand(n,n));
// [u,s,v] = svd(a);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// test_val = t1all;
//@}
//@{ test_svd9.m
//% Test the compact version of the svd with float vectors
//function test_val = test_svd9
// t1all = 1;
// p = [];
// for n=2:4:100
// a = float(rand(n,2*n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// p(n) = err;
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = float(rand(2*n,n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = float(rand(n,n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// test_val = t1all;
//@}
//@{ test_svd10.m
//% Test the compact version of the svd with complex vectors
//function test_val = test_svd10
// t1all = 1;
// p = [];
// for n=2:4:100
// a = float(rand(n,2*n))+i*float(rand(n,2*n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// p(n) = err;
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = float(rand(2*n,n))+i*float(rand(2*n,n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:4:100
// a = float(rand(n,n))+i*float(rand(n,n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*feps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// test_val = t1all;
//@}
//@{ test_svd11.m
//% Test the compact version of the svd with double vectors
//function test_val = test_svd11
// t1all = 1;
// p = [];
// for n=2:100
// a = double(rand(n,2*n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// p(n) = err;
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:100
// a = double(rand(2*n,n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:100
// a = double(rand(n,n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// test_val = t1all;
//@}
//@{ test_svd12.m
//% Test the compact version of the svd with dcomplex vectors
//function test_val = test_svd12
// t1all = 1;
// p = [];
// for n=2:100
// a = double(rand(n,2*n))+i*double(rand(n,2*n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// p(n) = err;
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:100
// a = double(rand(2*n,n))+i*double(rand(2*n,n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// for n=2:100
// a = double(rand(n,n))+i*double(rand(n,n));
// [u,s,v] = svd(a,0);
// emat = abs(a-u*s*v');
// err = max(emat(:));
// bnd = 10*max(abs(diag(s)))*eps*n;
// t1 = (err < bnd);
// if (~t1) printf('test failed: er = %e bnd = %e (num %d)\n',err,bnd,n); end;
// t1all = t1all & t1;
// end
// test_val = t1all;
//@}
//!
ArrayVector SVDFunction(int nargout, const ArrayVector& arg) {
if (arg.size() > 2)
throw Exception("svd function takes at most two arguments");
if (arg.size() < 1)
throw Exception("svd function requries at least one argument - the matrix to decompose");
Array A(arg[0]);
bool compactform = false;
if (arg.size() > 1) {
Array flag(arg[1]);
if (flag.getContentsAsIntegerScalar() == 0)
compactform = true;
}
// Test for numeric
if (A.isReferenceType())
throw Exception("Cannot apply svd to reference types.");
if (!A.is2D())
throw Exception("Cannot apply matrix operations to N-Dimensional arrays.");
if (A.anyNotFinite())
throw Exception("SVD only defined for matrices with finite entries.");
int nrows = A.getDimensionLength(0);
int ncols = A.getDimensionLength(1);
Class Aclass(A.dataClass());
if (Aclass < FM_FLOAT) {
A.promoteType(FM_DOUBLE);
Aclass = FM_DOUBLE;
}
bool computevectors;
computevectors = (nargout>1);
ArrayVector retval;
switch (Aclass) {
case FM_FLOAT:
{
int mindim;
mindim = (nrows < ncols) ? nrows : ncols;
// A temporary vector to store the singular values
float *svals = (float*) Malloc(mindim*sizeof(float));
// A temporary vector to store the left singular vectors
float *umat = NULL;
// A temporary vector to store the right singular vectors
float *vtmat = NULL;
if (computevectors)
if (!compactform) {
umat = (float*) Malloc(nrows*nrows*sizeof(float));
vtmat = (float*) Malloc(ncols*ncols*sizeof(float));
} else {
umat = (float*) Malloc(nrows*mindim*sizeof(float));
vtmat = (float*) Malloc(ncols*mindim*sizeof(float));
}
floatSVD(nrows,ncols,umat,vtmat,svals,
(float*) A.getReadWriteDataPointer(),
compactform, computevectors);
// Always transfer the singular values into an Array
Dimensions dim;
if (!computevectors) {
dim.set(0,mindim); dim.set(1,1);
retval.push_back(Array(FM_FLOAT,dim,svals));
} else {
if (!compactform) {
dim.set(0,nrows); dim.set(1,nrows);
retval.push_back(Array(FM_FLOAT,dim,umat));
dim.set(0,nrows); dim.set(1,ncols);
float *smat = (float*) Malloc(nrows*ncols*sizeof(float));
for (int i=0;i<mindim;i++)
smat[i+i*nrows] = svals[i];
retval.push_back(Array(FM_FLOAT,dim,smat));
dim.set(0,ncols); dim.set(1,ncols);
Array Utrans(FM_FLOAT,dim,vtmat);
Utrans.transpose();
retval.push_back(Utrans);
Free(svals);
} else {
dim.set(0,nrows); dim.set(1,mindim);
retval.push_back(Array(FM_FLOAT,dim,umat));
dim.set(0,mindim); dim.set(1,mindim);
float *smat = (float*) Malloc(mindim*mindim*sizeof(float));
for (int i=0;i<mindim;i++)
smat[i+i*mindim] = svals[i];
retval.push_back(Array(FM_FLOAT,dim,smat));
dim.set(0,mindim); dim.set(1,ncols);
Array Utrans(FM_FLOAT,dim,vtmat);
Utrans.transpose();
retval.push_back(Utrans);
Free(svals);
}
}
break;
}
case FM_DOUBLE:
{
int mindim;
mindim = (nrows < ncols) ? nrows : ncols;
// A temporary vector to store the singular values
double *svals = (double*) Malloc(mindim*sizeof(double));
// A temporary vector to store the left singular vectors
double *umat = NULL;
// A temporary vector to store the right singular vectors
double *vtmat = NULL;
if (computevectors)
if (!compactform) {
umat = (double*) Malloc(nrows*nrows*sizeof(double));
vtmat = (double*) Malloc(ncols*ncols*sizeof(double));
} else {
umat = (double*) Malloc(nrows*mindim*sizeof(double));
vtmat = (double*) Malloc(ncols*mindim*sizeof(double));
}
doubleSVD(nrows,ncols,umat,vtmat,svals,
(double*) A.getReadWriteDataPointer(),
compactform, computevectors);
// Always transfer the singular values into an Array
Dimensions dim;
if (!computevectors) {
dim.set(0,mindim); dim.set(1,1);
retval.push_back(Array(FM_DOUBLE,dim,svals));
Free(umat);
Free(vtmat);
} else {
if (!compactform) {
dim.set(0,nrows); dim.set(1,nrows);
retval.push_back(Array(FM_DOUBLE,dim,umat));
dim.set(0,nrows); dim.set(1,ncols);
double *smat = (double*) Malloc(nrows*ncols*sizeof(double));
for (int i=0;i<mindim;i++)
smat[i+i*nrows] = svals[i];
retval.push_back(Array(FM_DOUBLE,dim,smat));
dim.set(0,ncols); dim.set(1,ncols);
Array Utrans(FM_DOUBLE,dim,vtmat);
Utrans.transpose();
retval.push_back(Utrans);
Free(svals);
} else {
dim.set(0,nrows); dim.set(1,mindim);
retval.push_back(Array(FM_DOUBLE,dim,umat));
dim.set(0,mindim); dim.set(1,mindim);
double *smat = (double*) Malloc(mindim*mindim*sizeof(double));
for (int i=0;i<mindim;i++)
smat[i+i*mindim] = svals[i];
retval.push_back(Array(FM_DOUBLE,dim,smat));
dim.set(0,mindim); dim.set(1,ncols);
Array Utrans(FM_DOUBLE,dim,vtmat);
Utrans.transpose();
retval.push_back(Utrans);
Free(svals);
}
}
break;
}
case FM_COMPLEX:
{
int mindim;
mindim = (nrows < ncols) ? nrows : ncols;
// A temporary vector to store the singular values
float *svals = (float*) Malloc(mindim*sizeof(float));
// A temporary vector to store the left singular vectors
float *umat = NULL;
// A temporary vector to store the right singular vectors
float *vtmat = NULL;
if (computevectors)
if (!compactform) {
umat = (float*) Malloc(2*nrows*nrows*sizeof(float));
vtmat = (float*) Malloc(2*ncols*ncols*sizeof(float));
} else {
umat = (float*) Malloc(2*nrows*mindim*sizeof(float));
vtmat = (float*) Malloc(2*ncols*mindim*sizeof(float));
}
complexSVD(nrows,ncols,umat,vtmat,svals,
(float*) A.getReadWriteDataPointer(),
compactform, computevectors);
// Always transfer the singular values into an Array
Dimensions dim;
if (!computevectors) {
dim.set(0,mindim); dim.set(1,1);
retval.push_back(Array(FM_FLOAT,dim,svals));
Free(umat);
Free(vtmat);
} else {
if (!compactform) {
dim.set(0,nrows); dim.set(1,nrows);
retval.push_back(Array(FM_COMPLEX,dim,umat));
dim.set(0,nrows); dim.set(1,ncols);
float *smat = (float*) Malloc(nrows*ncols*sizeof(float));
for (int i=0;i<mindim;i++)
smat[i+i*nrows] = svals[i];
retval.push_back(Array(FM_FLOAT,dim,smat));
dim.set(0,ncols); dim.set(1,ncols);
Array Utrans(FM_COMPLEX,dim,vtmat);
Utrans.hermitian();
retval.push_back(Utrans);
Free(svals);
} else {
dim.set(0,nrows); dim.set(1,mindim);
retval.push_back(Array(FM_COMPLEX,dim,umat));
dim.set(0,mindim); dim.set(1,mindim);
float *smat = (float*) Malloc(mindim*mindim*sizeof(float));
for (int i=0;i<mindim;i++)
smat[i+i*mindim] = svals[i];
retval.push_back(Array(FM_FLOAT,dim,smat));
dim.set(0,mindim); dim.set(1,ncols);
Array Utrans(FM_COMPLEX,dim,vtmat);
Utrans.hermitian();
retval.push_back(Utrans);
Free(svals);
}
}
break;
}
case FM_DCOMPLEX:
{
int mindim;
mindim = (nrows < ncols) ? nrows : ncols;
// A temporary vector to store the singular values
double *svals = (double*) Malloc(mindim*sizeof(double));
// A temporary vector to store the left singular vectors
double *umat = NULL;
// A temporary vector to store the right singular vectors
double *vtmat = NULL;
if (computevectors)
if (!compactform) {
umat = (double*) Malloc(2*nrows*nrows*sizeof(double));
vtmat = (double*) Malloc(2*ncols*ncols*sizeof(double));
} else {
umat = (double*) Malloc(2*nrows*mindim*sizeof(double));
vtmat = (double*) Malloc(2*ncols*mindim*sizeof(double));
}
dcomplexSVD(nrows,ncols,umat,vtmat,svals,
(double*) A.getReadWriteDataPointer(),
compactform, computevectors);
// Always transfer the singular values into an Array
Dimensions dim;
if (!computevectors) {
dim.set(0,mindim); dim.set(1,1);
retval.push_back(Array(FM_DOUBLE,dim,svals));
Free(umat);
Free(vtmat);
} else {
if (!compactform) {
dim.set(0,nrows); dim.set(1,nrows);
retval.push_back(Array(FM_DCOMPLEX,dim,umat));
dim.set(0,nrows); dim.set(1,ncols);
double *smat = (double*) Malloc(nrows*ncols*sizeof(double));
for (int i=0;i<mindim;i++)
smat[i+i*nrows] = svals[i];
retval.push_back(Array(FM_DOUBLE,dim,smat));
dim.set(0,ncols); dim.set(1,ncols);
Array Utrans(FM_DCOMPLEX,dim,vtmat);
Utrans.hermitian();
retval.push_back(Utrans);
Free(svals);
} else {
dim.set(0,nrows); dim.set(1,mindim);
retval.push_back(Array(FM_DCOMPLEX,dim,umat));
dim.set(0,mindim); dim.set(1,mindim);
double *smat = (double*) Malloc(mindim*mindim*sizeof(double));
for (int i=0;i<mindim;i++)
smat[i+i*mindim] = svals[i];
retval.push_back(Array(FM_DOUBLE,dim,smat));
dim.set(0,mindim); dim.set(1,ncols);
Array Utrans(FM_DCOMPLEX,dim,vtmat);
Utrans.hermitian();
retval.push_back(Utrans);
Free(svals);
}
}
}
break;
}
return retval;
}
//!
//@Module LASTERR Retrieve Last Error Message
//@@Section FLOW
//@@Usage
//Either returns or sets the last error message. The
//general syntax for its use is either
//@[
// msg = lasterr
//@]
//which returns the last error message that occured, or
//@[
// lasterr(msg)
//@]
//which sets the contents of the last error message.
//@@Example
//Here is an example of using the @|error| function to
//set the last error, and then retrieving it using
//lasterr.
//@<
//try; error('Test error message'); catch; end;
//lasterr
//@>
//Or equivalently, using the second form:
//@<
//lasterr('Test message');
//lasterr
//@>
//!
ArrayVector LasterrFunction(int nargout, const ArrayVector& arg,
Interpreter* eval) {
ArrayVector retval;
if (arg.size() == 0) {
Array A = Array::stringConstructor(eval->getLastErrorString());
retval.push_back(A);
} else {
Array tmp(arg[0]);
eval->setLastErrorString(tmp.getContentsAsString());
}
return retval;
}
//!
//@Module SIMKEYS Simulate Keypresses from the User
//@@Section FREEMAT
//@@Usage
//This routine simulates keystrokes from the user on FreeMat.
//The general syntax for its use is
//@[
// otext = simkeys(text)
//@]
//where @|text| is a string to simulate as input to the console.
//The output of the commands are captured and returned in the
//string @|otext|. This is primarily used by the testing
//infrastructure.
//!
ArrayVector SimKeysFunction(int nargout, const ArrayVector& arg,
Interpreter* eval) {
if (arg.size() == 0)
throw Exception("simkeys requires at least one argument (the cell array of strings to simulate)");
eval->clearCaptureString();
eval->setCaptureState(true);
if (arg[0].dataClass() != FM_CELL_ARRAY)
throw Exception("simkeys requires a cell array of strings");
const Array *dp = (const Array *) arg[0].getDataPointer();
for (int i=0;i<arg[0].getLength();i++)
eval->ExecuteLine(ArrayToString(dp[i]));
eval->ExecuteLine("quit\n");
try {
while(1)
eval->evalCLI();
} catch (InterpreterContinueException& e) {
} catch (InterpreterBreakException& e) {
} catch (InterpreterReturnException& e) {
} catch (InterpreterRetallException& e) {
} catch (InterpreterQuitException& e) {
}
eval->setCaptureState(false);
return ArrayVector() << Array::stringConstructor(eval->getCaptureString());
}
//!
//@Module ERRORCOUNT Retrieve the Error Counter for the Interpreter
//@@Section FREEMAT
//@@Usage
//This routine retrieves the internal counter for the interpreter,
//and resets it to zero. The general syntax for its use is
//@[
// count = errorcount
//@]
//!
ArrayVector ErrorCountFunction(int nargout, const ArrayVector& arg,
Interpreter* eval) {
return ArrayVector() << Array::int32Constructor(eval->getErrorCount());
}
//!
//@Module DIAG Diagonal Matrix Construction/Extraction
//@@Section ARRAY
//@@Usage
//The @|diag| function is used to either construct a
//diagonal matrix from a vector, or return the diagonal
//elements of a matrix as a vector. The general syntax
//for its use is
//@[
// y = diag(x,n)
//@]
//If @|x| is a matrix, then @|y| returns the @|n|-th
//diagonal. If @|n| is omitted, it is assumed to be
//zero. Conversely, if @|x| is a vector, then @|y|
//is a matrix with @|x| set to the @|n|-th diagonal.
//@@Examples
//Here is an example of @|diag| being used to extract
//a diagonal from a matrix.
//@<
//A = int32(10*rand(4,5))
//diag(A)
//diag(A,1)
//@>
//Here is an example of the second form of @|diag|, being
//used to construct a diagonal matrix.
//@<
//x = int32(10*rand(1,3))
//diag(x)
//diag(x,-1)
//@>
//@@Tests
//@{ test_diag1.m
//% Test the diagonal extraction function
//function test_val = test_diag1
//a = [1,2,3,4;5,6,7,8;9,10,11,12];
//b = diag(a);
//test_val = test(b == [1;6;11]);
//@}
//@{ test_diag2.m
//% Test the diagonal extraction function with a non-zero diagonal
//function test_val = test_diag2
//a = [1,2,3,4;5,6,7,8;9,10,11,12];
//b = diag(a,1);
//test_val = test(b == [2;7;12]);
//@}
//@{ test_diag3.m
//% Test the diagonal creation function
//function test_val = test_diag3
//a = [2,3];
//b = diag(a);
//test_val = test(b == [2,0;0,3]);
//@}
//@{ test_diag4.m
//% Test the diagonal creation function with a non-zero diagonal
//function test_val = test_diag4
//a = [2,3];
//b = diag(a,-1);
//test_val = test(b == [0,0,0;2,0,0;0,3,0]);
//@}
//@{ test_diag5.m
//% Test the diagonal creation function with no arguments (bug 1620051)
//function test_val = test_diag5
//test_val = 1;
//try
// b = diag;
//catch
// test_val = 1;
//end
//@}
//@@Tests
//@{ test_sparse74.m
//% Test sparse matrix array diagonal extraction
//function x = test_sparse74
//[yi1,zi1] = sparse_test_mat('int32',300,400);
//[yf1,zf1] = sparse_test_mat('float',300,400);
//[yd1,zd1] = sparse_test_mat('double',300,400);
//[yc1,zc1] = sparse_test_mat('complex',300,400);
//[yz1,zz1] = sparse_test_mat('dcomplex',300,400);
//x = testeq(diag(yi1,30),diag(zi1,30)) & testeq(diag(yf1,30),diag(zf1,30)) & testeq(diag(yd1,30),diag(zd1,30)) & testeq(diag(yc1,30),diag(zc1,30)) & testeq(diag(yz1,30),diag(zz1,30));
//@}
//!
ArrayVector DiagFunction(int nargout, const ArrayVector& arg) {
Array a;
Array b;
Array c;
ArrayVector retval;
int32 *dp;
int diagonalOrder;
// First, get the diagonal order, and synthesize it if it was
// not supplied
if (arg.size() == 0)
throw Exception("diag requires at least one argument.\n");
if (arg.size() == 1)
diagonalOrder = 0;
else {
b = arg[1];
if (!b.isScalar())
throw Exception("second argument must be a scalar.\n");
b.promoteType(FM_INT32);
dp = (int32 *) b.getDataPointer();
diagonalOrder = dp[0];
}
// Next, make sure the first argument is 2D
a = arg[0];
if (!a.is2D())
throw Exception("First argument to 'diag' function must be 2D.\n");
// Case 1 - if the number of columns in a is 1, then this is a diagonal
// constructor call.
if ((a.getDimensionLength(1) == 1) || (a.getDimensionLength(0) == 1))
c = Array::diagonalConstructor(a,diagonalOrder);
else
c = a.getDiagonal(diagonalOrder);
retval.push_back(c);
return retval;
}
//!
//@Module ISEMPTY Test For Variable Empty
//@@Section INSPECTION
//@@Usage
//The @|isempty| function returns a boolean that indicates
//if the argument variable is empty or not. The general
//syntax for its use is
//@[
// y = isempty(x).
//@]
//@@Examples
//Here are some examples of the @|isempty| function
//@<
//a = []
//isempty(a)
//b = 1:3
//isempty(b)
//@>
//Note that if the variable is not defined, @|isempty|
//does not return true.
//@<1
//clear x
//isempty(x)
//@>
//@@Tests
//@{ test_empty.m
//% Test the arithmetic operators with empty arguments
//function test_val = test_empty
//test_val = isempty([]+[]);
//test_val = test_val & isempty([]-[]);
//test_val = test_val & isempty(-[]);
//test_val = test_val & isempty([]*[]);
//test_val = test_val & isempty([]/[]);
//test_val = test_val & isempty([]\[]);
//test_val = test_val & isempty([].*[]);
//test_val = test_val & isempty([]./[]);
//test_val = test_val & isempty([].\[]);
//test_val = test_val & isempty([]^[]);
//test_val = test_val & isempty([].^[]);
//test_val = test_val & isempty([]>[]);
//test_val = test_val & isempty([]>=[]);
//test_val = test_val & isempty([]<[]);
//test_val = test_val & isempty([]<=[]);
//@}
//!
ArrayVector IsEmptyFunction(int nargout, const ArrayVector& arg) {
if (arg.size() < 1)
throw Exception("isempty function requires at least input argument");
ArrayVector retval;
retval.push_back(Array::logicalConstructor(arg[0].isEmpty()));
return retval;
}
//!
//@Module SPONES Sparse Ones Function
//@@Section SPARSE
//@@Usage
//Returns a sparse @|float| matrix with ones where the argument
//matrix has nonzero values. The general syntax for it is
//@[
// y = spones(x)
//@]
//where @|x| is a matrix (it may be full or sparse). The output
//matrix @|y| is the same size as @|x|, has type @|float|, and contains
//ones in the nonzero positions of @|x|.
//@@Examples
//Here are some examples of the @|spones| function
//@<
//a = [1,0,3,0,5;0,0,2,3,0;1,0,0,0,1]
//b = spones(a)
//full(b)
//@>
//!
ArrayVector SponesFunction(int nargout, const ArrayVector& arg) {
if (arg.size() < 1)
throw Exception("spones function requires a sparse matrix template argument");
Array tmp(arg[0]);
if (tmp.isEmpty())
return singleArrayVector(Array::emptyConstructor());
if(tmp.isReferenceType())
throw Exception("spones function requires a numeric sparse matrix argument");
tmp.makeSparse();
if (!tmp.sparse())
throw Exception("spones function requires a sparse matrix template argument");
return singleArrayVector(Array::Array(FM_FLOAT,Dimensions(tmp.getDimensionLength(0),tmp.getDimensionLength(1)),SparseOnesFunc(tmp.dataClass(),tmp.getDimensionLength(0),tmp.getDimensionLength(1),tmp.getSparseDataPointer()),true));
}
//!
//@Module WARNING Emits a Warning Message
//@@Section FLOW
//@@Usage
//The @|warning| function causes a warning message to be
//sent to the user. The general syntax for its use is
//@[
// warning(s)
//@]
//where @|s| is the string message containing the warning.
//!
ArrayVector WarningFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
if (arg.size() == 0)
throw Exception("Not enough inputs to warning function");
if (!(arg[0].isString()))
throw Exception("Input to error function must be a string");
eval->warningMessage(arg[0].getContentsAsString());
return ArrayVector();
}
//!
//@Module VERSION The Current Version Number
//@@Section FREEMAT
//@@Usage
//The @|version| function returns the current version number for
//FreeMat (as a string). The general syntax for its use is
//@[
// v = version
//@]
//@@Example
//The current version of FreeMat is
//@<
//version
//@>
//!
ArrayVector VersionFunction(int nargout, const ArrayVector& arg) {
return ArrayVector() << Array::stringConstructor(VERSION);
}
//!
//@Module VERSTRING The Current Version String
//@@Section FREEMAT
//@@Usage
//The @|verstring| function returns the current version string for
//FreeMat. The general syntax for its use is
//@[
// version = verstring
//@]
//@@Example
//The current version of FreeMat is
//@<
//verstring
//@>
//!
ArrayVector VerStringFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
return ArrayVector() << Array::stringConstructor(Interpreter::getVersionString());
}
//!
//@Module ERROR Causes an Error Condition Raised
//@@Section FLOW
//@@Usage
//The @|error| function causes an error condition (exception
//to be raised). The general syntax for its use is
//@[
// error(s),
//@]
//where @|s| is the string message describing the error. The
//@|error| function is usually used in conjunction with @|try|
//and @|catch| to provide error handling. If the string @|s|,
//then (to conform to the MATLAB API), @|error| does nothing.
//@@Example
//Here is a simple example of an @|error| being issued by a function
//@|evenoddtest|:
//@{ evenoddtest.m
//function evenoddtest(n)
// if (n==0)
// error('zero is neither even nor odd');
// elseif (~isa(n,'int32'))
// error('expecting integer argument');
// end;
// if (n==int32(n/2)*2)
// printf('%d is even\n',n);
// else
// printf('%d is odd\n',n);
// end
//@}
//The normal command line prompt @|-->| simply prints the error
//that occured.
//@<2
//evenoddtest(4)
//evenoddtest(5)
//evenoddtest(0)
//evenoddtest(pi)
//@>
//!
ArrayVector ErrorFunction(int nargout, const ArrayVector& arg) {
if (arg.size() == 0)
return ArrayVector();
string etxt(ArrayToString(arg[0]));
if (!etxt.empty())
throw Exception(etxt);
return ArrayVector();
}
//!
//@Module EVAL Evaluate a String
//@@Section FREEMAT
//@@Usage
//The @|eval| function evaluates a string. The general syntax
//for its use is
//@[
// eval(s)
//@]
//where @|s| is the string to evaluate. If @|s| is an expression
//(instead of a set of statements), you can assign the output
//of the @|eval| call to one or more variables, via
//@[
// x = eval(s)
// [x,y,z] = eval(s)
//@]
//
//Another form of @|eval| allows you to specify an expression or
//set of statements to execute if an error occurs. In this
//form, the syntax for @|eval| is
//@[
// eval(try_clause,catch_clause),
//@]
//or with return values,
//@[
// x = eval(try_clause,catch_clause)
// [x,y,z] = eval(try_clause,catch_clause)
//@]
//These later forms are useful for specifying defaults. Note that
//both the @|try_clause| and @|catch_clause| must be expressions,
//as the equivalent code is
//@[
// try
// [x,y,z] = try_clause
// catch
// [x,y,z] = catch_clause
// end
//@]
//so that the assignment must make sense in both cases.
//@@Example
//Here are some examples of @|eval| being used.
//@<
//eval('a = 32')
//b = eval('a')
//@>
//The primary use of the @|eval| statement is to enable construction
//of expressions at run time.
//@<
//s = ['b = a' ' + 2']
//eval(s)
//@>
//Here we demonstrate the use of the catch-clause to provide a
//default value
//@<
//a = 32
//b = eval('a','1')
//b = eval('z','a+1')
//@>
//Note that in the second case, @|b| takes the value of 33, indicating
//that the evaluation of the first expression failed (because @|z| is
//not defined).
//!
static string PrePendCallVars(string line, int nargout) {
char buf[4096];
char *gp = buf;
if (nargout > 1)
*gp++ = '[';
for (int i=0;i<nargout-1;i++) {
sprintf(gp,"t___%d,",i);
gp += strlen(gp);
}
sprintf(gp,"t___%d",nargout-1);
gp += strlen(gp);
if (nargout > 1)
sprintf(gp,"] = %s;\n",line.c_str());
else
sprintf(gp," = %s;\n",line.c_str());
return buf;
}
//
//static - having a static here seems to break cross compiling
//
ArrayVector RetrieveCallVars(Interpreter *eval, int nargout) {
ArrayVector retval;
for (int i=0;i<nargout;i++) {
char tname[4096];
Array tval;
sprintf(tname,"t___%d",i);
ArrayReference ptr = eval->getContext()->lookupVariable(tname);
if (!ptr.valid())
tval = Array::emptyConstructor();
else
tval = *ptr;
eval->getContext()->deleteVariable(tname);
retval.push_back(tval);
}
return retval;
}
ArrayVector EvalTryFunction(int nargout, const ArrayVector& arg, Interpreter* eval, int popSpec) {
if (nargout > 0) {
string try_line = arg[0].getContentsAsString();
string try_buf = PrePendCallVars(try_line,nargout);
string catch_line = arg[1].getContentsAsString();
string catch_buf = PrePendCallVars(catch_line,nargout);
ArrayVector retval;
bool autostop;
autostop = eval->AutoStop();
eval->AutoStop(false);
try {
eval->getContext()->bypassScope(popSpec);
eval->evaluateString(try_buf,true);
retval = RetrieveCallVars(eval,nargout);
eval->getContext()->restoreBypassedScopes();
} catch (Exception &e) {
eval->getContext()->restoreBypassedScopes();
eval->evaluateString(catch_buf,false);
retval = RetrieveCallVars(eval,nargout);
}
eval->AutoStop(autostop);
return retval;
} else {
string try_line = arg[0].getContentsAsString();
string catch_line = arg[1].getContentsAsString();
string try_buf = try_line + "\n";
string catch_buf = catch_line + "\n";
bool autostop;
autostop = eval->AutoStop();
eval->AutoStop(false);
try {
eval->getContext()->bypassScope(popSpec);
eval->evaluateString(try_buf,true);
eval->getContext()->restoreBypassedScopes();
} catch (Exception &e) {
eval->getContext()->restoreBypassedScopes();
eval->evaluateString(catch_buf,false);
}
eval->AutoStop(autostop);
return ArrayVector();
}
}
ArrayVector EvalNoTryFunction(int nargout, const ArrayVector& arg, Interpreter* eval, int popSpec) {
if (nargout > 0) {
string line = arg[0].getContentsAsString();
string buf = PrePendCallVars(line,nargout);
eval->getContext()->bypassScope(popSpec);
eval->evaluateString(buf);
ArrayVector retval(RetrieveCallVars(eval,nargout));
eval->getContext()->restoreBypassedScopes();
return retval;
} else {
string line = arg[0].getContentsAsString();
string buf = line + "\n";
eval->getContext()->bypassScope(popSpec);
eval->evaluateString(buf);
eval->getContext()->restoreBypassedScopes();
return ArrayVector();
}
}
ArrayVector EvalFunction(int nargout, const ArrayVector& arg,Interpreter* eval){
if (arg.size() == 0)
throw Exception("eval function takes at least one argument - the string to execute");
if (arg.size() == 2)
return EvalTryFunction(nargout, arg, eval, 0);
return EvalNoTryFunction(nargout, arg, eval, 0);
}
//!
//@Module DIARY Create a Log File of Console
//@@Section FREEMAT
//@@Usage
//The @|diary| function controls the creation of a log file that duplicates
//the text that would normally appear on the console.
//The simplest syntax for the command is simply:
//@[
// diary
//@]
//which toggles the current state of the diary command. You can also explicitly
//set the state of the diary command via the syntax
//@[
// diary off
//@]
//or
//@[
// diary on
//@]
//To specify a filename for the log (other than the default of @|diary|), you
//can use the form:
//@[
// diary filename
//@]
//or
//@[
// diary('filename')
//@]
//which activates the diary with an output filename of @|filename|. Note that the
//@|diary| command is thread specific, but that the output is appended to a given
//file. That means that if you call @|diary| with the same filename on multiple
//threads, their outputs will be intermingled in the log file (just as on the console).
//Because the @|diary| state is tied to individual threads, you cannot retrieve the
//current diary state using the @|get(0,'Diary')| syntax from MATLAB. Instead, you
//must call the @|diary| function with no inputs and one output:
//@[
// state = diary
//@]
//which returns a logical @|1| if the output of the current thread is currently going to
//a diary, and a logical @|0| if not.
//!
ArrayVector DiaryFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
if (nargout == 1) {
if (arg.size() > 0)
throw Exception("diary function with an assigned return value (i.e, 'x=diary') does not support any arguments");
return ArrayVector() << Array::logicalConstructor(eval->getDiaryState());
}
if (arg.size() == 0) {
eval->setDiaryState(!eval->getDiaryState());
return ArrayVector();
}
string diaryString(ArrayToString(arg[0]));
if (diaryString == "on")
eval->setDiaryState(true);
else if (diaryString == "off")
eval->setDiaryState(false);
else {
eval->setDiaryFilename(diaryString);
eval->setDiaryState(true);
}
return ArrayVector();
}
//!
//@Module EVALIN Evaluate a String in Workspace
//@@Section FREEMAT
//@@Usage
//The @|evalin| function is similar to the @|eval| function, with an additional
//argument up front that indicates the workspace that the expressions are to
//be evaluated in. The various syntaxes for @|evalin| are:
//@[
// evalin(workspace,expression)
// x = evalin(workspace,expression)
// [x,y,z] = evalin(workspace,expression)
// evalin(workspace,try_clause,catch_clause)
// x = evalin(workspace,try_clause,catch_clause)
// [x,y,z] = evalin(workspace,try_clause,catch_clause)
//@]
//The argument @|workspace| must be either 'caller' or 'base'. If it is
//'caller', then the expression is evaluated in the caller's work space.
//That does not mean the caller of @|evalin|, but the caller of the current
//function or script. On the other hand if the argument is 'base', then
//the expression is evaluated in the base work space. See @|eval| for
//details on the use of each variation.
//!
ArrayVector EvalInFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
if (arg.size() < 2)
throw Exception("evalin function requires a workspace (scope) specifier (either 'caller' or 'base') and an expression to evaluate");
Array spec(arg[0]);
string spec_str = spec.getContentsAsString();
int popspec = 0;
if (spec_str=="base")
popspec = -1;
else if (spec_str=="caller")
popspec = 1;
else
throw Exception("evalin function requires the first argument to be either 'caller' or 'base'");
ArrayVector argcopy(arg);
argcopy.pop_front();
if (arg.size() == 3)
return EvalTryFunction(nargout,argcopy,eval,popspec);
else
return EvalNoTryFunction(nargout,argcopy,eval,popspec);
}
//!
//@Module ASSIGNIN Assign Variable in Workspace
//@@Section FREEMAT
//@@Usage
//The @|assignin| function allows you to assign a value to a variable
//in either the callers work space or the base work space. The syntax
//for @|assignin| is
//@[
// assignin(workspace,variablename,value)
//@]
//The argument @|workspace| must be either 'caller' or 'base'. If it is
//'caller' then the variable is assigned in the caller's work space.
//That does not mean the caller of @|assignin|, but the caller of the
//current function or script. On the other hand if the argument is 'base',
//then the assignment is done in the base work space. Note that the
//variable is created if it does not already exist.
//!
ArrayVector AssignInFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
if (arg.size() < 3)
throw Exception("assignin function requires a workspace (scope) specifier (either 'caller' or 'base') a variable name and a value to assign");
Array spec(arg[0]);
string spec_str = spec.getContentsAsString();
int popspec = 0;
if (spec_str=="base")
popspec = -1;
else if (spec_str=="caller")
popspec = 1;
else
throw Exception("evalin function requires the first argument to be either 'caller' or 'base'");
string varname = ArrayToString(arg[1]);
Array varvalue = arg[2];
eval->getContext()->bypassScope(popspec);
eval->getContext()->insertVariable(varname,varvalue);
eval->getContext()->restoreBypassedScopes();
return ArrayVector();
}
//!
//@Module QUIET Control the Verbosity of the Interpreter
//@@Section FREEMAT
//@@Usage
//The @|quiet| function controls how verbose the interpreter
//is when executing code. The syntax for the function is
//@[
// quiet flag
//@]
//where @|flag| is one of
//\begin{itemize}
//\item @|'normal'| - normal output from the interpreter
//\item @|'quiet'| - only intentional output (e.g. @|printf| calls and
//@|disp| calls) is printed. The output of expressions that are not
//terminated in semicolons are not printed.
//\item @|'silent'| - nothing is printed to the output.
//\end{itemize}
//The @|quiet| command also returns the current quiet flag.
//!
ArrayVector QuietFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
if (arg.size() > 0) {
string qtype(arg[0].getContentsAsStringUpper());
if (qtype == "NORMAL")
eval->setQuietLevel(0);
else if (qtype == "QUIET")
eval->setQuietLevel(1);
else if (qtype == "SILENT")
eval->setQuietLevel(2);
else
throw Exception("quiet function takes one argument - the quiet level (normal, quiet, or silent) as a string");
}
string rtype;
if (eval->getQuietLevel() == 0)
rtype = "normal";
else if (eval->getQuietLevel() == 1)
rtype = "quiet";
else if (eval->getQuietLevel() == 2)
rtype = "silent";
return ArrayVector() << Array::stringConstructor(rtype);
}
//!
//@Module SOURCE Execute an Arbitrary File
//@@Section FREEMAT
//@@Usage
//The @|source| function executes the contents of the given
//filename one line at a time (as if it had been typed at
//the @|-->| prompt). The @|source| function syntax is
//@[
// source(filename)
//@]
//where @|filename| is a @|string| containing the name of
//the file to process.
//@@Example
//First, we write some commands to a file (note that it does
//not end in the usual @|.m| extension):
//@{ source_test
//a = 32;
//b = a;
//printf('a is %d and b is %d\n',a,b);
//@}
//Now we source the resulting file.
//@<
//clear a b
//source source_test
//@>
//@@Tests
//@{ source_test_script.m
//n = 1;
//n = n + 1;
//@}
//@{ test_source.m
//% Check that the source function does not double-execute the last line of the script
//function test_val = test_source
//myloc = which('test_source');
//[path,name,sfx] = fileparts(myloc);
//source([path '/source_test_script.m'])
//test_val = test(n == 2);
//@}
//!
ArrayVector SourceFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
if (arg.size() != 1)
throw Exception("source function takes exactly one argument - the filename of the script to execute");
string filename = arg[0].getContentsAsString();
QFile fp(QString::fromStdString(filename));
if (!fp.open(QFile::ReadOnly))
throw Exception(std::string("unable to open file ") + filename + " for reading");
QTextStream fstr(&fp);
QString scriptText(fstr.readAll());
Scanner S(scriptText.toStdString(),filename);
Parser P(S);
tree pcode = P.Process();
if (pcode.is(TOK_FUNCTION_DEFS))
throw Exception("only scripts can be source-ed, not functions");
tree code = pcode.first();
eval->block(code);
return ArrayVector();
}
//!
//@Module DBDELETE Delete a Breakpoint
//@@Section DEBUG
//@@Usage
//The @|dbdelete| function deletes a breakpoint. The syntax
//for the @|dbdelete| function is
//@[
// dbdelete(num)
//@]
//where @|num| is the number of the breakpoint to delete.
//!
ArrayVector DbDeleteFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
if (arg.size() < 1)
throw Exception("dbdelete requires an argument (number of breakpoint to delete)");
Array tmp(arg[0]);
int bpnum;
bpnum = tmp.getContentsAsIntegerScalar();
eval->deleteBreakpoint(bpnum);
return ArrayVector();
}
//!
//@Module DBLIST List Breakpoints
//@@Section DEBUG
//@@Usage
//List the current set of breakpoints. The syntax for the
//@|dblist| is simply
//@[
// dblist
//@]
//!
ArrayVector DbListFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
eval->listBreakpoints();
return ArrayVector();
}
//!
//@Module DBSTEP Step N Statements
//@@Section DEBUG
//@@Usage
//Step @|N| statements during debug mode. The synax for this is
//either
//@[
// dbstep(N)
//@]
//to step @|N| statements, or
//@[
// dbstep
//@]
//to step one statement.
//!
// ArrayVector DbStepFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
// int linesToSkip;
// if (arg.size() == 0)
// linesToSkip = 1;
// else {
// Array tmp(arg[0]);
// linesToSkip = tmp.getContentsAsIntegerScalar();
// }
// eval->dbstep(linesToSkip);
// return ArrayVector();
// }
//!
//@Module DBSTOP
//@@Section DEBUG
//@@Usage
//Set a breakpoint. The syntax for this is:
//@[
// dbstop(funcname,linenumber)
//@]
//where @|funcname| is the name of the function where we want
//to set the breakpoint, and @|linenumber| is the line number.
//!
ArrayVector DbStopFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
if (arg.size() < 2)
throw Exception("dbstop function requires at least two arguments");
if (!(arg[0].isString()))
throw Exception("first argument to dbstop must be the name of routine where to stop");
string cname = ArrayToString(arg[0]);
int line = ArrayToInt32(arg[1]);
eval->addBreakpoint(cname,line);
return ArrayVector();
}
ArrayVector FdumpFunction(int nargout, const ArrayVector& arg,Interpreter* eval){
if (arg.size() == 0)
throw Exception("fdump function requires at least one argument");
if (!(arg[0].isString()))
throw Exception("first argument to fdump must be the name of a function (i.e., a string)");
string fname = arg[0].getContentsAsString();
Context *context = eval->getContext();
FuncPtr funcDef;
if (!context->lookupFunction(fname,funcDef))
throw Exception(std::string("function ") + fname + " undefined!");
funcDef->updateCode(eval);
funcDef->printMe(eval);
return ArrayVector();
}
//!
//@Module BUILTIN Evaulate Builtin Function
//@@Section FREEMAT
//@@Usage
//The @|builtin| function evaluates a built in function
//with the given name, bypassing any overloaded functions.
//The syntax of @|builtin| is
//@[
// [y1,y2,...,yn] = builtin(fname,x1,x2,...,xm)
//@]
//where @|fname| is the name of the function to call. Apart
//from the fact that @|fname| must be a string, and that @|builtin|
//always calls the non-overloaded method, it operates exactly like
//@|feval|. Note that unlike MATLAB, @|builtin| does not force
//evaluation to an actual compiled function. It simply subverts
//the activation of overloaded method calls.
//!
ArrayVector BuiltinFunction(int nargout, const ArrayVector& arg,Interpreter* eval){
if (arg.size() == 0)
throw Exception("builtin function requires at least one argument");
if (!(arg[0].isString()))
throw Exception("first argument to builtin must be the name of a function (i.e., a string)");
FuncPtr funcDef;
string fname = arg[0].getContentsAsString();
Context *context = eval->getContext();
if (!context->lookupFunction(fname,funcDef))
throw Exception(std::string("function ") + fname + " undefined!");
funcDef->updateCode(eval);
if (funcDef->scriptFlag)
throw Exception("cannot use feval on a script");
ArrayVector newarg(arg);
newarg.pop_front();
bool flagsave = eval->getStopOverload();
eval->setStopOverload(true);
ArrayVector tmp(funcDef->evaluateFunction(eval,newarg,nargout));
eval->setStopOverload(flagsave);
return tmp;
}
//!
//@Module FEVAL Evaluate a Function
//@@Section FREEMAT
//@@Usage
//The @|feval| function executes a function using its name.
//The syntax of @|feval| is
//@[
// [y1,y2,...,yn] = feval(f,x1,x2,...,xm)
//@]
//where @|f| is the name of the function to evaluate, and
//@|xi| are the arguments to the function, and @|yi| are the
//return values.
//
//Alternately, @|f| can be a function handle to a function
//(see the section on @|function handles| for more information).
//
//Finally, FreeMat also supports @|f| being a user defined class
//in which case it will atttempt to invoke the @|subsref| method
//of the class.
//@@Example
//Here is an example of using @|feval| to call the @|cos|
//function indirectly.
//@<
//feval('cos',pi/4)
//@>
//Now, we call it through a function handle
//@<
//c = @cos
//feval(c,pi/4)
//@>
//Here we construct an inline object (which is a user-defined class)
//and use @|feval| to call it
//@<
//afunc = inline('cos(t)+sin(t)','t')
//feval(afunc,pi)
//afunc(pi)
//@>
//In both cases, (the @|feval| call and the direct invokation), FreeMat
//calls the @|subsref| method of the class, which computes the requested
//function.
//@@Tests
//@$"y = feval(@cos,pi)","-1","close"
//@$"y = feval('cos',pi)","-1","close"
//@$"y = feval(inline('cos(t)'),pi)","-1","close"
//!
ArrayVector FevalFunction(int nargout, const ArrayVector& arg,Interpreter* eval){
if (arg.size() == 0)
throw Exception("feval function requires at least one argument");
if (!(arg[0].isString()) && (arg[0].dataClass() != FM_FUNCPTR_ARRAY))
throw Exception("first argument to feval must be the name of a function (i.e., a string) a function handle, or a user defined class");
FuncPtr funcDef;
if (arg[0].isString()) {
string fname = arg[0].getContentsAsString();
Context *context = eval->getContext();
if (!context->lookupFunction(fname,funcDef))
throw Exception(std::string("function ") + fname + " undefined!");
} else {
if (!arg[0].isScalar())
throw Exception("function handle argument to feval must be a scalar");
const FuncPtr *fp = (const FuncPtr *) arg[0].getDataPointer();
funcDef = (FuncPtr) fp[0];
if (!funcDef) return ArrayVector();
}
funcDef->updateCode(eval);
if (funcDef->scriptFlag)
throw Exception("cannot use feval on a script");
ArrayVector newarg(arg);
newarg.pop_front();
return(funcDef->evaluateFunction(eval,newarg,nargout));
}
//!
//@Module DOCLI Start a Command Line Interface
//@@Section FREEMAT
//@@Usage
//The @|docli| function is the main function that you interact with
//when you run FreeMat. I am not sure why you would want to use
//it, but hey - its there if you want to use it.
//!
//!
//@Module STARTUP Startup Script
//@@Section FREEMAT
//@@Usage
//Upon starting, FreeMat searches for a script names @|startup.m|, and
//if it finds it, it executes it. This script can be in the current
//directory, or on the FreeMat path (set using @|setpath|). The contents
//of startup.m must be a valid script (not a function).
//!
ArrayVector DoCLIFunction(int nargout, const ArrayVector& arg, Interpreter* eval) {
Context *context = eval->getContext();
FuncPtr funcDef;
if (context->lookupFunction("startup",funcDef)) {
funcDef->updateCode(eval);
if (funcDef->scriptFlag) {
try {
eval->block(((MFunctionDef*)funcDef)->code);
} catch (Exception& e) {
eval->errorMessage("Startup script error:\n" + e.getMessageCopy());
}
} else {
eval->outputMessage(string("startup.m must be a script"));
}
}
eval->doCLI();
return ArrayVector();
}
//!
//@Module REPMAT Array Replication Function
//@@Section ARRAY
//@@Usage
//The @|repmat| function replicates an array the specified
//number of times. The source and destination arrays may
//be multidimensional. There are three distinct syntaxes for
//the @|repmap| function. The first form:
//@[
// y = repmat(x,n)
//@]
//replicates the array @|x| on an @|n-times-n| tiling, to create
//a matrix @|y| that has @|n| times as many rows and columns
//as @|x|. The output @|y| will match @|x| in all remaining
//dimensions. The second form is
//@[
// y = repmat(x,m,n)
//@]
//And creates a tiling of @|x| with @|m| copies of @|x| in the
//row direction, and @|n| copies of @|x| in the column direction.
//The final form is the most general
//@[
// y = repmat(x,[m n p...])
//@]
//where the supplied vector indicates the replication factor in
//each dimension.
//@@Example
//Here is an example of using the @|repmat| function to replicate
//a row 5 times. Note that the same effect can be accomplished
//(although somewhat less efficiently) by a multiplication.
//@<
//x = [1 2 3 4]
//y = repmat(x,[5,1])
//@>
//The @|repmat| function can also be used to create a matrix of scalars
//or to provide replication in arbitrary dimensions. Here we use it to
//replicate a 2D matrix into a 3D volume.
//@<
//x = [1 2;3 4]
//y = repmat(x,[1,1,3])
//@>
//!
#define MAX(a,b) ((a) > (b) ? (a) : (b))
ArrayVector RepMatFunction(int nargout, const ArrayVector& arg) {
int i, j, k;
if (arg.size() < 2)
throw Exception("repmat function requires at least two arguments");
Array x(arg[0]);
if (x.isEmpty()) return singleArrayVector(Array::emptyConstructor());
Dimensions repcount;
// Case 1, look for a scalar second argument
if ((arg.size() == 2) && (arg[1].isScalar())) {
Array t(arg[1]);
repcount.set(0,t.getContentsAsIntegerScalar());
repcount.set(1,t.getContentsAsIntegerScalar());
}
// Case 2, look for two scalar arguments
else if ((arg.size() == 3) && (arg[1].isScalar()) && (arg[2].isScalar())) {
Array t(arg[1]);
repcount.set(0,t.getContentsAsIntegerScalar());
t = arg[2];
repcount.set(1,t.getContentsAsIntegerScalar());
}
// Case 3, look for a vector second argument
else {
if (arg.size() > 2) throw Exception("unrecognized form of arguments for repmat function");
Array t(arg[1]);
t.promoteType(FM_INT32);
if (t.getLength() > maxDims) throw Exception("replication request exceeds maxDims constant - either rebuild FreeMat with a higher maxDims constant, or shorten the requested replication array");
int32 *dp = (int32*) t.getDataPointer();
for (i=0;i<t.getLength();i++)
repcount.set(i,dp[i]);
}
// Check for a valid replication count
for (i=0;i<repcount.getLength();i++)
if (repcount.get(i) < 0) throw Exception("negative replication counts not allowed in argument to repmat function");
// All is peachy. Allocate an output array of sufficient size.
Dimensions originalSize(x.dimensions());
Dimensions outdims;
int outdim;
outdim = MAX(repcount.getLength(),originalSize.getLength());
for (i=0;i<outdim;i++)
outdims.set(i,originalSize.get(i)*repcount.get(i));
outdims.simplify();
void *dp = Array::allocateArray(x.dataClass(),
outdims.getElementCount(),
x.fieldNames());
// Copy can work by pushing or by pulling. I have opted for
// pushing, because we can push a column at a time, which might
// be slightly more efficient.
int colsize = originalSize.get(0);
int outcolsize = outdims.get(0);
int colcount = originalSize.getElementCount()/colsize;
// copySelect stores which copy we are pushing.
Dimensions copySelect(outdim);
// anchor is used to calculate where this copy lands in the output matrix
Dimensions anchor(outdim);
// sourceAddress is used to track which column we are pushing in the
// source matrix
Dimensions sourceAddress(outdims.getLength());
int destanchor;
int copyCount;
copyCount = repcount.getElementCount();
for (i=0;i<copyCount;i++) {
// Reset the source address
sourceAddress.zeroOut();
// Next, we loop over the columns of the source matrix
for (j=0;j<colcount;j++) {
// We can calculate the anchor of this copy by multiplying the source
// address by the copySelect vector
for (k=0;k<outdim;k++)
anchor.set(k,copySelect.get(k)*originalSize.get(k)+sourceAddress.get(k));
// Now, we map this to a point in the destination array
destanchor = outdims.mapPoint(anchor);
// And copy the elements
x.copyElements(j*colsize,dp,destanchor,colsize);
// Now increment the source address
sourceAddress.incrementModulo(originalSize,1);
}
copySelect.incrementModulo(repcount,0);
}
ArrayVector retval;
retval.push_back(Array(x.dataClass(),outdims,dp,false,x.fieldNames()));
return retval;
}
//!
//@Module SYSTEM Call an External Program
//@@Section OS
//@@Usage
//The @|system| function allows you to call an external
//program from within FreeMat, and capture the output.
//The syntax of the @|system| function is
//@[
// y = system(cmd)
//@]
//where @|cmd| is the command to execute. The return
//array @|y| is of type @|cell-array|, where each entry
//in the array corresponds to a line from the output.
//@@Example
//Here is an example of calling the @|ls| function (the
//list files function under Un*x-like operating system).
//@<
//y = system('ls m*.m')
//y{1}
//@>
//!
ArrayVector SystemFunction(int nargout, const ArrayVector& arg) {
if (arg.size() != 1)
throw Exception("System function takes one string argument");
string systemArg = arg[0].getContentsAsString();
ArrayVector retval;
if (systemArg.size() == 0)
return retval;
stringVector cp(DoSystemCallCaptured(systemArg));
Array* np = new Array[cp.size()];
for (int i=0;i<cp.size();i++)
np[i] = Array::stringConstructor(cp[i]);
Dimensions dim(cp.size(),1);
Array ret(FM_CELL_ARRAY,dim,np);
retval.push_back(ret);
return retval;
}
//!
//@Module PERMUTE Array Permutation Function
//@@Section ARRAY
//@@Usage
//The @|permute| function rearranges the contents of an array according
//to the specified permutation vector. The syntax for its use is
//@[
// y = permute(x,p)
//@]
//where @|p| is a permutation vector - i.e., a vector containing the
//integers @|1...ndims(x)| each occuring exactly once. The resulting
//array @|y| contains the same data as the array @|x|, but ordered
//according to the permutation. This function is a generalization of
//the matrix transpose operation.
//@@Example
//Here we use @|permute| to transpose a simple matrix (note that permute
//also works for sparse matrices):
//@<
//A = [1,2;4,5]
//permute(A,[2,1])
//A'
//@>
//Now we permute a larger n-dimensional array:
//@<
//A = randn(13,5,7,2);
//size(A)
//B = permute(A,[3,4,2,1]);
//size(B)
//@>
//@@Tests
//@$"y=permute([1,2;3,4],[2,1])","[1,3;2,4]","exact"
//@$"y=size(permute(randn(13,5,7,2),[3,4,2,1]))","[7,2,5,13]","exact"
//@{ test_permute1.m
//function test_val = test_permute1
//z = rand(3,5,2,4,7);
//perm = [3,5,1,4,2];
//sizez = size(z);
//y = permute(z,perm);
//sizey = size(y);
//test_val = all(sizey == sizez(perm));
//@}
//!
ArrayVector PermuteFunction(int nargout, const ArrayVector& arg) {
if (arg.size() < 2) throw Exception("permute requires 2 inputs, the array to permute, and the permutation vector");
Array A(arg[0]);
Array permutation(arg[1]);
permutation.promoteType(FM_INT32);
int Adims = A.dimensions().getLength();
if (permutation.getLength() != A.dimensions().getLength())
throw Exception("permutation vector must contain as many elements as the array to permute has dimensions");
// Check that it is, in fact a permutation
MemBlock<bool> p(Adims);
bool *d = &p;
for (int i=0;i<Adims;i++) d[i] = false;
const int32* dp = (const int32*) permutation.getDataPointer();
for (int i=0;i<Adims;i++) {
if ((dp[i] < 1) || (dp[i] > Adims))
throw Exception("permutation vector elements are limited to 1..ndims(A), where A is the array to permute");
d[dp[i]-1] = true;
}
// Check that all are covered
for (int i=0;i<Adims;i++)
if (!d[i]) throw Exception("second argument to permute function is not a permutation (no duplicates allowed)");
A.permute(dp);
return ArrayVector() << A;
}
//!
//@Module CONV2 Matrix Convolution
//@@Section SIGNAL
//@@Usage
//The @|conv2| function performs a two-dimensional convolution of
//matrix arguments. The syntax for its use is
//@[
// Z = conv2(X,Y)
//@]
//which performs the full 2-D convolution of @|X| and @|Y|. If the
//input matrices are of size @|[xm,xn]| and @|[ym,yn]| respectively,
//then the output is of size @|[xm+ym-1,xn+yn-1]|. Another form is
//@[
// Z = conv2(hcol,hrow,X)
//@]
//where @|hcol| and @|hrow| are vectors. In this form, @|conv2|
//first convolves @|Y| along the columns with @|hcol|, and then
//convolves @|Y| along the rows with @|hrow|. This is equivalent
//to @|conv2(hcol(:)*hrow(:)',Y)|.
//
//You can also provide an optional @|shape| argument to @|conv2|
//via either
//@[
// Z = conv2(X,Y,'shape')
// Z = conv2(hcol,hrow,X,'shape')
//@]
//where @|shape| is one of the following strings
//\begin{itemize}
//\item @|'full'| - compute the full convolution result - this is the default if no @|shape| argument is provided.
//\item @|'same'| - returns the central part of the result that is the same size as @|X|.
//\item @|'valid'| - returns the portion of the convolution that is computed without the zero-padded edges. In this situation, @|Z| has
//size @|[xm-ym+1,xn-yn+1]| when @|xm>=ym| and @|xn>=yn|. Otherwise
//@|conv2| returns an empty matrix.
//\end{itemize}
//@@Function Internals
//The convolution is computed explicitly using the definition:
//\[
// Z(m,n) = \sum_{k} \sum_{j} X(k,j) Y(m-k,n-j)
//\]
//If the full output is requested, then @|m| ranges over @|0 <= m < xm+ym-1|
//and @|n| ranges over @|0 <= n < xn+yn-1|. For the case where @|shape|
//is @|'same'|, the output ranges over @|(ym-1)/2 <= m < xm + (ym-1)/2|
//and @|(yn-1)/2 <= n < xn + (yn-1)/2|.
//!
// Summation limits...
// sum a(i,j) b(m-i,n-j)
// Need
// 0 <= m-i <= Bm-1
// Can rewrite this as:
// -m <= -i <= Bm-1-m
// m >= i >= m+1-Bm
// m+1-Bm <= i <= m
// 0 <= n-j <= Bn-1
// n >= j >= n+1-Bn
// n+1-Bn <= j <= n
// And
// 0 <= i <= Am-1
// 0 <= j <= An-1
// So
// m+1-Bm <= i <= m
// 0 <= i <= Am-1
// max(0,m+1-Bm) <= i <= min(m,Am-1)
//
// n+1-Bn <= j <= n
// 0 <= j <= An-1
// max(0,n+1-Bn) <= j <= min(n,An-1)
//
//
// max(0,m+1-Bm) <= i <= min(m,Am-1)
// max(0,n+1-Bn) <= j <= min(n,An-1)
template <class T>
static void Conv2MainReal(T* C, const T* A, const T*B,
int Am, int An, int Bm, int Bn,
int Cm, int Cn, int Cm_offset, int Cn_offset) {
for (int n=0;n<Cn;n++)
for (int m=0;m<Cm;m++) {
T accum = 0;
int iMin, iMax;
int jMin, jMax;
iMin = std::max(0,m+Cm_offset-Bm+1);
iMax = std::min(Am-1,m+Cm_offset);
jMin = std::max(0,n+Cn_offset-Bn+1);
jMax = std::min(An-1,n+Cn_offset);
for (int j=jMin;j<=jMax;j++)
for (int i=iMin;i<=iMax;i++)
accum += A[i+j*Am]*B[(m+Cm_offset-i)+(n+Cn_offset-j)*Bm];
C[m+n*Cm] = accum;
}
}
template <class T>
static void Conv2MainComplex(T* C, const T* A, const T*B,
int Am, int An, int Bm, int Bn,
int Cm, int Cn, int Cm_offset, int Cn_offset) {
for (int n=0;n<Cn;n++)
for (int m=0;m<Cm;m++) {
T accum_r = 0;
T accum_i = 0;
int iMin, iMax;
int jMin, jMax;
iMin = std::max(0,m+Cm_offset-Bm+1);
iMax = std::min(Am-1,m+Cm_offset);
jMin = std::max(0,n+Cn_offset-Bn+1);
jMax = std::min(An-1,n+Cn_offset);
for (int j=jMin;j<=jMax;j++)
for (int i=iMin;i<=iMax;i++) {
accum_r += A[2*(i+j*Am)]*
B[2*((m+Cm_offset-i)+(n+Cn_offset-j)*Bm)]
- A[2*(i+j*Am)+1]*B[2*((m+Cm_offset-i)+(n+Cn_offset-j)*Bm)+1];
accum_i += A[2*(i+j*Am)]*
B[2*((m+Cm_offset-i)+(n+Cn_offset-j)*Bm)+1]
- A[2*(i+j*Am)+1]*B[2*((m+Cm_offset-i)+(n+Cn_offset-j)*Bm)];
}
C[2*(m+n*Cm)] = accum_r;
C[2*(m+n*Cm)] = accum_i;
}
}
Array Conv2FunctionDispatch(Array X,Array Y,int Cm,int Cn,
int Cm_offset, int Cn_offset) {
switch (X.dataClass()) {
case FM_FLOAT: {
float *cp = (float*) Array::allocateArray(FM_FLOAT,Cm*Cn);
Conv2MainReal<float>(cp,
(const float*) X.getDataPointer(),
(const float*) Y.getDataPointer(),
X.getDimensionLength(0),
X.getDimensionLength(1),
Y.getDimensionLength(0),
Y.getDimensionLength(1),
Cm, Cn, Cm_offset, Cn_offset);
return Array(FM_FLOAT,Dimensions(Cm,Cn),cp);
}
case FM_DOUBLE: {
double *cp = (double*) Array::allocateArray(FM_DOUBLE,Cm*Cn);
Conv2MainReal<double>(cp,
(const double*) X.getDataPointer(),
(const double*) Y.getDataPointer(),
X.getDimensionLength(0),
X.getDimensionLength(1),
Y.getDimensionLength(0),
Y.getDimensionLength(1),
Cm, Cn, Cm_offset, Cn_offset);
return Array(FM_DOUBLE,Dimensions(Cm,Cn),cp);
}
case FM_INT32: {
int32 *cp = (int32*) Array::allocateArray(FM_INT32,Cm*Cn);
Conv2MainReal<int32>(cp,
(const int32*) X.getDataPointer(),
(const int32*) Y.getDataPointer(),
X.getDimensionLength(0),
X.getDimensionLength(1),
Y.getDimensionLength(0),
Y.getDimensionLength(1),
Cm, Cn, Cm_offset, Cn_offset);
return Array(FM_INT32,Dimensions(Cm,Cn),cp);
}
case FM_INT64: {
int64 *cp = (int64*) Array::allocateArray(FM_INT64,Cm*Cn);
Conv2MainReal<int64>(cp,
(const int64*) X.getDataPointer(),
(const int64*) Y.getDataPointer(),
X.getDimensionLength(0),
X.getDimensionLength(1),
Y.getDimensionLength(0),
Y.getDimensionLength(1),
Cm, Cn, Cm_offset, Cn_offset);
return Array(FM_INT64,Dimensions(Cm,Cn),cp);
}
case FM_COMPLEX: {
float *cp = (float*) Array::allocateArray(FM_COMPLEX,Cm*Cn);
Conv2MainComplex<float>(cp,
(const float*) X.getDataPointer(),
(const float*) Y.getDataPointer(),
X.getDimensionLength(0),
X.getDimensionLength(1),
Y.getDimensionLength(0),
Y.getDimensionLength(1),
Cm, Cn, Cm_offset, Cn_offset);
return Array(FM_COMPLEX,Dimensions(Cm,Cn),cp);
}
case FM_DCOMPLEX: {
double *cp = (double*) Array::allocateArray(FM_DCOMPLEX,Cm*Cn);
Conv2MainComplex<double>(cp,
(const double*) X.getDataPointer(),
(const double*) Y.getDataPointer(),
X.getDimensionLength(0),
X.getDimensionLength(1),
Y.getDimensionLength(0),
Y.getDimensionLength(1),
Cm, Cn, Cm_offset, Cn_offset);
return Array(FM_DCOMPLEX,Dimensions(Cm,Cn),cp);
}
}
}
static ArrayVector Conv2FunctionFullXY(Array X, Array Y) {
int Cm, Cn, Cm_offset, Cn_offset;
Cm = X.getDimensionLength(0) + Y.getDimensionLength(0) - 1;
Cn = X.getDimensionLength(1) + Y.getDimensionLength(1) - 1;
Cm_offset = 0;
Cn_offset = 0;
return singleArrayVector(Conv2FunctionDispatch(X,Y,Cm,Cn,Cm_offset,Cn_offset));
}
static ArrayVector Conv2FunctionSameXY(Array X, Array Y) {
int Cm, Cn, Cm_offset, Cn_offset;
Cm = X.getDimensionLength(0);
Cn = X.getDimensionLength(1);
Cm_offset = (int) floor((double)((Y.getDimensionLength(0)-1)/2));
Cn_offset = (int) floor((double)((Y.getDimensionLength(1)-1)/2));
return singleArrayVector(Conv2FunctionDispatch(X,Y,Cm,Cn,Cm_offset,Cn_offset));
}
static ArrayVector Conv2FunctionValidXY(Array X, Array Y) {
int Cm, Cn, Cm_offset, Cn_offset;
Cm = X.getDimensionLength(0)-Y.getDimensionLength(0)+1;
Cn = X.getDimensionLength(1)-Y.getDimensionLength(1)+1;
if ((Cm <= 0) || (Cn <= 0))
return singleArrayVector(Array::emptyConstructor());
Cm_offset = Y.getDimensionLength(0)-1;
Cn_offset = Y.getDimensionLength(1)-1;
return singleArrayVector(Conv2FunctionDispatch(X,Y,Cm,Cn,Cm_offset,Cn_offset));
}
ArrayVector Conv2FunctionXY(Array X, Array Y, string type) {
// Check the arguments
if (X.isReferenceType() || Y.isReferenceType())
throw Exception("cannot apply conv2 to reference types.");
if (!X.is2D() || !Y.is2D())
throw Exception("arguments must be matrices, not n-dimensional arrays.");
TypeCheck(X,Y,false);
if (type == "FULL")
return Conv2FunctionFullXY(X,Y);
if (type == "SAME")
return Conv2FunctionSameXY(X,Y);
if (type == "VALID")
return Conv2FunctionValidXY(X,Y);
throw Exception("could not recognize the arguments to conv2");
}
ArrayVector Conv2FunctionRCX(Array hcol, Array hrow, Array X, string type) {
if (hcol.isReferenceType() || hrow.isReferenceType() ||
X.isReferenceType())
throw Exception("cannot apply conv2 to reference types.");
if (!hcol.is2D() || !hrow.is2D() || !X.is2D())
throw Exception("arguments must be matrices, not n-dimensional arrays.");
Dimensions t(hcol.getLength(),1);
hcol.reshape(t);
t = Dimensions(1,hrow.getLength());
hrow.reshape(t);
ArrayVector rvec;
rvec = Conv2FunctionXY(X,hcol,type);
rvec = Conv2FunctionXY(rvec.back(),hrow,type);
return rvec;
}
ArrayVector Conv2Function(int nargout, const ArrayVector& arg) {
// Call the right function based on the arguments
if (arg.size() < 2)
throw Exception("conv2 requires at least 2 arguments");
if (arg.size() == 2)
return Conv2FunctionXY(arg[0],arg[1],"FULL");
if ((arg.size() == 3) && (arg[2].isString()))
return Conv2FunctionXY(arg[0],arg[1],arg[2].getContentsAsStringUpper());
if (arg.size() == 3)
return Conv2FunctionRCX(arg[0],arg[1],arg[2],"FULL");
if ((arg.size() == 4) && (arg[3].isString()))
return Conv2FunctionRCX(arg[0],arg[1],arg[2],
arg[3].getContentsAsStringUpper());
throw Exception("could not recognize which form of conv2 was requested - check help conv2 for details");
}
syntax highlighted by Code2HTML, v. 0.9.1