文件列表:

linalg\LinAlg\LinAlg.h (57685, 1998-12-23)
linalg\LinAlg\LAStreams.h (27761, 1998-12-07)
linalg\LinAlg\Makefile (3506, 1998-12-28)
linalg\LinAlg\ali.cc (8581, 1998-12-02)
linalg\LinAlg\builtin.h (3013, 1998-12-21)
linalg\LinAlg\determinant.cc (6482, 1998-10-12)
linalg\LinAlg\fminbr.cc (5763, 1997-06-12)
linalg\LinAlg\hjmin.cc (6395, 1998-12-02)
linalg\LinAlg\math_num.h (3169, 1998-10-05)
linalg\LinAlg\matrix1.cc (20498, 1998-12-19)
linalg\LinAlg\matrix2.cc (5396, 1998-12-19)
linalg\LinAlg\matrix_inv.cc (5350, 1998-12-19)
linalg\LinAlg\matrix_sub.cc (7720, 1998-12-02)
linalg\LinAlg\minmax.h (2452, 1996-12-21)
linalg\LinAlg\myenv.cc (7364, 1998-12-28)
linalg\LinAlg\myenv.h (3393, 1998-12-21)
linalg\LinAlg\sample_adv.cc (3603, 1998-12-21)
linalg\LinAlg\sample_adv.dat (345, 1998-01-23)
linalg\LinAlg\std.h (1245, 1998-04-13)
linalg\LinAlg\svd.cc (24895, 1998-12-19)
linalg\LinAlg\svd.h (3602, 1998-12-19)
linalg\LinAlg\vali.cc (5181, 1998-01-03)
linalg\LinAlg\vali.dat (3588, 1998-12-26)
linalg\LinAlg\vector.cc (4206, 1997-12-30)
linalg\LinAlg\vfminbr.cc (2615, 1998-12-17)
linalg\LinAlg\vfminbr.dat (1125, 1998-12-26)
linalg\LinAlg\vhjmin.cc (6167, 1998-12-02)
linalg\LinAlg\vhjmin.dat (2057, 1998-12-26)
linalg\LinAlg\vmatrix.cc (10106, 1998-11-25)
linalg\LinAlg\vmatrix.dat (2838, 1998-12-26)
linalg\LinAlg\vmatrix1.cc (8540, 1998-10-26)
linalg\LinAlg\vmatrix1.dat (1192, 1998-12-26)
linalg\LinAlg\vmatrix2.cc (11132, 1998-12-21)
linalg\LinAlg\vmatrix2.dat (3390, 1998-12-26)
linalg\LinAlg\vslesing.cc (1976, 1998-12-14)
linalg\LinAlg\vslesing.dat (3526, 1998-12-26)
linalg\LinAlg\vsvd.cc (4508, 1998-12-14)
linalg\LinAlg\vsvd.dat (34097, 1998-12-26)
linalg\LinAlg\vvector.cc (19799, 1998-12-02)
... ...

Numerical Math Class Library version 4.3 ***** Version history is at the very end ***** Comments, questions, problem reports, etc. are all very welcome. Please mail me at oleg@pobox.com or oleg@acm.org or oleg@computer.org ***** Platforms I have personally compiled and tested LinAlg 4.3 on Linux 2.0.33 i586, gcc 2.7.2.1 HP 9000/770, HP-UX 10.10, gcc 2.8.1 (with optimization off) UltraSparcII, Solaris 2.6, gcc 2.8.1 (with optimization off) PowerMac 8500/132, BeOS Release 4, Metrowerk's CodeWarrior C++ Pentium, WinNT 4.0, Visual C++ 5.0 The previous (4.1) version also works on SunSparc20/Solaris 2.4, gcc 2.7.2 and SunWorkshopPro SGI, gcc 2.7.2 and SGI's native compiler PowerMac 7100/80, 8500/132, Metrowerk's CodeWarrior C++, v. 11-12 DEC Alpha (Digital UNIX Version 5.60), gcc 2.7.2 Note: if g++ version 2.8.1 is used to compile LinAlg, optimization must be turned off. Otherwise, the compiler crashes with internal compiler errors. Other compilers as well as gcc v2.7.2 don't seem to have this problem. ***** Verification files: vmatrix vvector vmatrix1 vmatrix2 vlastreams vali vhjmin vfminbr vzeroin vsvd vslesing Don't forget to compile and run them, see comments in the Makefile for details. The verification code checks to see that all the functions in this package have compiled and run well. The validation code can also serve as an illustration of how package's classes and functions may be employed. ***** Highlights and idioms 1. Data classes and view classes Objects of a class Matrix are the only ones that _own_ data, that is, 2D grids of REALs. All other classes provide different views of/into the grids: in 2D, in 1D (a view of all matrix elements linearly arranged in a column major order), a view of a single matrix row, column, or the diagonal, a view of a rectangular subgrid, etc. The views are designed to be as lightweight as possible: most of the view functionality could be inlined. That is why views do not have any virtual functions. View classes are supposed to be final, in Java parlance: It makes little sense to inherit from them. Indeed, views are designed to provide general and many special ways of accessing matrix elements, as safely and efficiently as possible. Views turn out to be so flexible that many former Matrix methods can now be implemented outside of the Matrix class. These methods can do their job without a privileged access to Matrix' internals while not sacrificing performance. Examples of such methods: comparison of two matrices, multiplication of a matrix by a diagonal of another matrix. 2. Matrix streams Matrix streams provide a sequential view/access to a matrix or its parts. Many of Linear Algebra algorithms actually require only sequential access to a matrix or its rows/columns, which is simpler and faster than a random access. That's why LinAlg stresses streams. There are two kinds of sequential views: ElementWise streams are transient views that exist only for the duration of an elementwise action. Hence these views don't require a "locking" of a matrix. On the other hand, LAStream and the ilk are more permanent views. An application traverses the stream at its own convenience, bookmarking certain spots and possibly returning to them later. An LAStream does assert a shared lock on the matrix, to prevent the grid of matrix values from moving or disappearing. Again, LAStreamIn etc. are meant to denote a stream under user's control, which the user can create and play with for some time, getting elements, moving the current position, etc. That is, LAStreamIn and the Matrix object from which the stream was created (whose view the stream is) can coexist. This arrangement is dangerous as one can potentially dispose of the container before all of its views are closed. To guard against this, a view asserts a "shared lock" on a matrix grid. Any operation that disposes of the grid or potentially moves it in memory (like resize_to()) checks this lock. Actually, it tries to assert an exclusive lock. Exclusive locking of a matrix object will always fail if there are any outstanding shared locks. In contrast, ElementWiseConst and ElementWise are meant to be transient, created only for the duration of a requested operation. That's why it is syntactically impossible to dispose of the matrix that is being viewed through a ElementWiseConst or ElementWise object. These objects do not have any public constructors. This prevents a user from building the view on his own (and forgetting to delete it). Therefore creation of ElementWiseConst objects does not require asserting any locks, which makes this class rather "light". Matrix streams may stride a matrix by an arbitrary amount. This allows traversing of a matrix along the diagonal, by columns, by rows, etc. Streams can be constructed of a Matrix itself, or from other matrix views (MatrixColumn, MatrixRow, MatrixDiag). In the latter case, the streams are confined only to specific portions of the matrix. Many methods and functions have been re-written to take advantage of the streams. Notable examples: // Vector norms are computed via the streams: inline double Vector::norm_1(void) const { return of_every(*this).sum_abs(); } inline double Vector::norm_2_sqr(void) const { return of_every(*this).sum_squares(); } inline double Vector::norm_inf(void) const { return of_every(*this).max_abs(); } The methods ElementWiseConst::sum_abs(), sum_squares(), and max_abs() can also be applied to two collections. In this case, they work through element-by-element differences between the collections. An example of using LAStreams: // Aitken-Lagrange interpolation (see ali.cc) double ALInterp::interpolate() { LAStreamIn args(arg); // arg and val are vectors LAStreamOut vals(val); register double valp = vals.peek(); // The best approximation found so far register double diffp = DBL_MAX; // abs(valp - prev. to valp) // Compute the j-th row of the Aitken scheme and // place it in the 'val' array for(register int j=2; j<=val.q_upb(); j++) { register double argj = (args.get(), args.peek()); register REAL& valj = (vals.get(), vals.peek()); args.rewind(); vals.rewind(); // rewind the vector streams for(register int i=1; i<=j-1; i++) { double argi = args.get(); valj = ( vals.get()*(q-argj) - valj*(q-argi) ) / (argi - argj); } .... } return valp; } Note, all the streams are light-weight: they use no heap storage and leave no garbage. All the operations in the snippets above are inlined. A stride of a stride stream doesn't have to be an even multiple of the number of elements in the corresponding collection, as the following not-so-trivial example shows. It places a 'value' at the _anti_-diagonal of a matrix m: m.clear(); LAStrideStreamOut m_str(m,m.q_nrows()-1); m_str.get(); // Skip the first element for(register int i=0; ij) inline void SVD::rotate(Matrix& U, const int i, const int j, const double cos_ph, const double sin_ph) { MatrixColumn Ui(U,i), Uj(U,j); for(register int l=1; l<=U.q_row_upb(); l++) { REAL& Uil = Ui(l); REAL& Ujl = Uj(l); const REAL Ujl_was = Ujl; Ujl = cos_ph*Ujl_was + sin_ph*Uil; Uil = -sin_ph*Ujl_was + cos_ph*Uil; } } The new version of the package, v4.3, rewrites this code in the following way, avoiding random access to Matrix elements (and incident range checks for MatrixColumns' indices). inline void SVD::rotate(Matrix& U, const int i, const int j, const double cos_ph, const double sin_ph) { LAStreamOut Ui(MatrixColumn (U,i)); LAStreamOut Uj(MatrixColumn (U,j)); while( !Ui.eof() ) { REAL& Uil = Ui.get(); REAL& Ujl = Uj.get(); const REAL Ujl_was = Ujl; Ujl = cos_ph*Ujl_was + sin_ph*Uil; Uil = -sin_ph*Ujl_was + cos_ph*Uil; } } LinAlg's implementation and validation code provides many more examples of flexibility, clarity, and efficiency of streams. 4. Subranging a stream LinAlg 4.3 implements subranging of a stream: creation of a new stream that spans over a part of another. The latter can be either a stride 1 or an arbitrary stride stream. The new stream may not _extend_ the old one: subranging is therefore a safe operation. A child stream is always confined within its father's boundaries. A substream, once created, lives independently of its parent: either of them can go out of scope without affecting the other. To create a substream, one has to specify its range as a [lwb,upb] pair: an all-inclusive range. The 'lwb' may also be given as -IRange::INF, and 'upb' as IRange::INF. These special values are interpreted intelligently. Given an output stream, one may create either an immutable, input substream, or allow modifications to a part of the original stream. Obviously, given a read-only stream, only read-only substreams may be created: the compiler will see to that at _compile_ time. File sample_ult.cc provides a good example of using substreams to efficiently reflect the upper triangle of a square matrix onto the lower one (yielding a symmetric matrix). See vlastreams.cc and SVD for more extensive examples of subranging. As a good illustration to the power and efficiency of streams, consider the following snippet from SVD::left_householder() (file svd.cc). In LinAlg 3.2, the snippet was programmed as: for(j=i+1; j<=N; j++) // Transform i+1:N columns of A { MatrixColumn Aj(A,j); REAL factor = 0; for(k=i; k<=M; k++) factor += UPi(k) * Aj(k); // Compute UPi' * A[,j] factor /= beta; for(k=i; k<=M; k++) Aj(k) -= UPi(k) * factor; } Version 4.3 uses substreams (which span only a part of matrix' columns): IRange range = IRange::from(i - A.q_col_lwb()); LAStreamOut UPi_str(MatrixColumn(A,i),range); ... for(register int j=i+1; j<=N; j++) // Transform i+1:N columns of A { LAStreamOut Aj_str(MatrixColumn(A,j),range); REAL factor = 0; while( !UPi_str.eof() ) factor += UPi_str.get() * Aj_str.get(); // Compute UPi' * A[,j] factor /= beta; for(UPi_str.rewind(), Aj_str.rewind(); !UPi_str.eof(); ) Aj_str.get() -= UPi_str.get() * factor; UPi_str.rewind(); } 5. Streams over an arbitrary rectangular block of a matrix LinAlg 4.3 permits LAstreams that span an arbitrary rectangular block of a matrix (including the whole matrix, a single matrix element, a matrix row or a column, or a part thereof). You simply specify a matrix, and the range of rows and columns to include in the stream. vlastreams.cc verification code shows very many examples. The following is a annotated snippet from vlastreams' output (vlastreams.dat) ---> Test LABlockStreams for 2:12x0:20 checking accessing of the whole matrix as a block stream... checking modifying the whole matrix as a block stream... # The first column of m checking LABlockStreams clipped as [-INF,INF] x [-INF,0] # The second column of m checking LABlockStreams clipped as [2,INF] x [1,1] # A part of the last column of m checking LABlockStreams clipped as [3,12] x [20,INF] # The first row of m checking LABlockStreams clipped as [-INF,2] x [0,INF] # The second row of m checking LABlockStreams clipped as [3,3] x [-INF,20] # A part of the last row of m checking LABlockStreams clipped as [12,INF] x [-INF,19] # The first matrix element checking LABlockStreams clipped as [2,2] x [0,0] # The last matrix element checking LABlockStreams clipped as [12,INF] x [20,INF] # A 2x3 block at the upper-right corner checking LABlockStreams clipped as [3,4] x [18,INF] # A 3x2 block at the lower-left corner checking LABlockStreams clipped as [9,11] x [2,3] checking modifying LABlockStreams clipped as [4,4] x [3,3] With block streams, it is trivial to assign a block of one matrix to a block of another matrix. See vlastreams.cc for examples. 6. Direct access to matrix elements In LinAlg 3.2 and earlier, a Matrix contained a column index. The index was used to speed up direct access to matrix elements (like in m(i,j)). The index was allocated on heap and initialized upon Matrix construction. In LinAlg 4+, a Matrix object no longer allocates any index. The Matrix class therefore does not allow direct access to matrix elements. This is done on purpose, to slim down matrices and make them easier to build. Many linear algebra operations actually require only sequential access to matrix elements. Thus the column index and other direct access facilities are often redundant. If one does need to access arbitrary elements of a matrix, the package offers several choices. One may use a special MatrixDA view, which is designed to provide an efficient direct access to a matrix grid. The MatrixDA view allocates and builds a column index, and uses it to efficiently compute a reference to the (i,j)-th element: Matrix ethc = m; MatrixDA eth(ethc); for(register int i=eth.q_row_lwb(); i<=eth.q_row_upb(); i++) for(register int j=eth.q_col_lwb(); j<=eth.q_col_upb(); j++) eth(i,j) *= v(j); see the validation code vmatrix*.cc, vvector.cc for more details. If creating of an index seems like overkill, one may use MatrixRow, matrixColumn or MatrixDiag. For example, the following statements are all equivalent: MatrixRow(m,2)(3) = M_PI; MatrixColumn(m,3)(2) = M_PI; MatrixDA ma(m); ma(2,3) = M_PI; (MatrixDA(m))(2,3) = M_PI; Unlike Matrix itself, a Vector and all matrix views (MatrixRow, MatrixColumn, MatrixDiag) allow direct access to their elements. In most of the cases, MatrixDA behaves like a matrix. On occasions when one needs to get hold of a real matrix, he can always use MatrixDA::ref() or MatrixDA::get_container() methods: Matrix mc(-1,10,1,20); MatrixDA m(mc); m.get_container().clear(); verify_element_value(m,0); Again, it is strongly encouraged to use stream-lined matrix operations as much as possible. As an example, finding of a matrix inverse and computing of a determinant are implemented in this package using only sequential matrix access. It is possible indeed to avoid direct access. 7. Const and Writable views There are two flavors of every matrix view: a const and a writable views, for example: ConstMatrixDA and MatrixDA, ConstMatrixColumn and MatrixColumn, ElementWiseConst and ElementWise, LAStreamIn and LAStreamOut. A const view can be constructed around a const matrix reference. A const view provides methods that do not alter the matrix: methods that query matrix elements, compare them, compute norms and other cumulative quantities (sum of squares, max absolute value, etc). A writable view is a subclass of the const view. That means that a writable view allows all the query/comparison operations that the corresponding const view implements. In addition, a writable view permits modification of matrix elements, via assignment or other related operations (+=, *=, etc). Needless to say one needs a non-const matrix reference to construct a writable view: cout << " compare (m = pattern^2)/pattern with pattern\n"; m = pattern; to_every(m1) = pattern; to_every(m).sqr(); assert( of_every(m).max_abs() == pattern*pattern ); to_every(m) /= of_every(m1); verify_element_value(m,pattern); to_every() makes a writable ElementWise view, while of_every() makes a ElementWiseConst view. 8. Never return non-trivial objects (matrices or vectors) Danger: For example, when the following snippet Matrix foo(const int n) { Matrix foom(n,n); fill_in(foom); return foom; } Matrix m = foo(5); runs, it constructs matrix foo:foom, copies it onto stack as a return value and destroys foo:foom. The return value (a matrix) from foo() is then copied over to m via a copy constructor. After that, the return value is destroyed. The matrix constructor is called 3 times and the destructor 2 times. For big matrices, the cost of multiple constructing/copying/destroying of objects may be very large. *Some* optimized compilers can do a little better. That still leaves at least two calls to the Matrix constructor. LazyMatrices (see below) can construct a Matrix m "inplace", with only a _single_ call to the constructor. 9. Lazy matrices Instead of returning an object return a "recipe" how to make it. The full matrix would be rolled out only when and where it is needed: Matrix haar = haar_matrix(5); haar_matrix() is a *class*, not a simple function. However similar this looks to returning of an object (see the note above), it is dramatically different. haar_matrix() constructs a LazyMatrix, an object just of a few bytes long. A special "Matrix(const LazyMatrix& recipe)" constructor follows the recipe and makes the matrix haar() right in place. No matrix element is moved whatsoever! Another example of matrix promises is REAL arr ... ...

近期下载者：

相关文件：

评论：[我要评论] [举报此文件]

收藏者：