From b2fc6c70434674d74551c3a6c01ffb3233499312 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 1 Jul 2013 22:34:11 +0000 Subject: Update version to 1.3 --- lib/include/rapidxml/rapidxml.hpp | 2596 +++++++++++++++++++++++++++ lib/include/rapidxml/rapidxml_print.hpp | 424 +++++ lib/include/tnt/jama_cholesky.h | 258 +++ lib/include/tnt/jama_eig.h | 1034 +++++++++++ lib/include/tnt/jama_lu.h | 323 ++++ lib/include/tnt/jama_qr.h | 326 ++++ lib/include/tnt/jama_svd.h | 543 ++++++ lib/include/tnt/tnt.h | 64 + lib/include/tnt/tnt_array1d.h | 278 +++ lib/include/tnt/tnt_array1d_utils.h | 230 +++ lib/include/tnt/tnt_array2d.h | 315 ++++ lib/include/tnt/tnt_array2d_utils.h | 287 +++ lib/include/tnt/tnt_array3d.h | 296 +++ lib/include/tnt/tnt_array3d_utils.h | 236 +++ lib/include/tnt/tnt_cmat.h | 580 ++++++ lib/include/tnt/tnt_fortran_array1d.h | 267 +++ lib/include/tnt/tnt_fortran_array1d_utils.h | 242 +++ lib/include/tnt/tnt_fortran_array2d.h | 225 +++ lib/include/tnt/tnt_fortran_array2d_utils.h | 236 +++ lib/include/tnt/tnt_fortran_array3d.h | 223 +++ lib/include/tnt/tnt_fortran_array3d_utils.h | 249 +++ lib/include/tnt/tnt_i_refvec.h | 243 +++ lib/include/tnt/tnt_math_utils.h | 34 + lib/include/tnt/tnt_sparse_matrix_csr.h | 103 ++ lib/include/tnt/tnt_stopwatch.h | 95 + lib/include/tnt/tnt_subscript.h | 54 + lib/include/tnt/tnt_vec.h | 404 +++++ lib/include/tnt/tnt_version.h | 39 + lib/licenses/rapidxml.txt | 52 + lib/licenses/tnt.txt | 14 + 30 files changed, 10270 insertions(+) create mode 100644 lib/include/rapidxml/rapidxml.hpp create mode 100644 lib/include/rapidxml/rapidxml_print.hpp create mode 100644 lib/include/tnt/jama_cholesky.h create mode 100644 lib/include/tnt/jama_eig.h create mode 100644 lib/include/tnt/jama_lu.h create mode 100644 lib/include/tnt/jama_qr.h create mode 100644 lib/include/tnt/jama_svd.h create mode 100644 lib/include/tnt/tnt.h create mode 100644 lib/include/tnt/tnt_array1d.h create mode 100644 lib/include/tnt/tnt_array1d_utils.h create mode 100644 lib/include/tnt/tnt_array2d.h create mode 100644 lib/include/tnt/tnt_array2d_utils.h create mode 100644 lib/include/tnt/tnt_array3d.h create mode 100644 lib/include/tnt/tnt_array3d_utils.h create mode 100644 lib/include/tnt/tnt_cmat.h create mode 100644 lib/include/tnt/tnt_fortran_array1d.h create mode 100644 lib/include/tnt/tnt_fortran_array1d_utils.h create mode 100644 lib/include/tnt/tnt_fortran_array2d.h create mode 100644 lib/include/tnt/tnt_fortran_array2d_utils.h create mode 100644 lib/include/tnt/tnt_fortran_array3d.h create mode 100644 lib/include/tnt/tnt_fortran_array3d_utils.h create mode 100644 lib/include/tnt/tnt_i_refvec.h create mode 100644 lib/include/tnt/tnt_math_utils.h create mode 100644 lib/include/tnt/tnt_sparse_matrix_csr.h create mode 100644 lib/include/tnt/tnt_stopwatch.h create mode 100644 lib/include/tnt/tnt_subscript.h create mode 100644 lib/include/tnt/tnt_vec.h create mode 100644 lib/include/tnt/tnt_version.h create mode 100644 lib/licenses/rapidxml.txt create mode 100644 lib/licenses/tnt.txt (limited to 'lib') diff --git a/lib/include/rapidxml/rapidxml.hpp b/lib/include/rapidxml/rapidxml.hpp new file mode 100644 index 0000000..ae91e08 --- /dev/null +++ b/lib/include/rapidxml/rapidxml.hpp @@ -0,0 +1,2596 @@ +#ifndef RAPIDXML_HPP_INCLUDED +#define RAPIDXML_HPP_INCLUDED + +// Copyright (C) 2006, 2009 Marcin Kalicinski +// Version 1.13 +// Revision $DateTime: 2009/05/13 01:46:17 $ +//! \file rapidxml.hpp This file contains rapidxml parser and DOM implementation + +// If standard library is disabled, user must provide implementations of required functions and typedefs +#if !defined(RAPIDXML_NO_STDLIB) + #include // For std::size_t + #include // For assert + #include // For placement new +#endif + +// On MSVC, disable "conditional expression is constant" warning (level 4). +// This warning is almost impossible to avoid with certain types of templated code +#ifdef _MSC_VER + #pragma warning(push) + #pragma warning(disable:4127) // Conditional expression is constant +#endif + +/////////////////////////////////////////////////////////////////////////// +// RAPIDXML_PARSE_ERROR + +#if defined(RAPIDXML_NO_EXCEPTIONS) + +#define RAPIDXML_PARSE_ERROR(what, where) { parse_error_handler(what, where); assert(0); } + +namespace rapidxml +{ + //! When exceptions are disabled by defining RAPIDXML_NO_EXCEPTIONS, + //! this function is called to notify user about the error. + //! It must be defined by the user. + //!

+ //! This function cannot return. If it does, the results are undefined. + //!

+ //! A very simple definition might look like that: + //! 
+    //! void %rapidxml::%parse_error_handler(const char *what, void *where)
+    //! {
+    //!     std::cout << "Parse error: " << what << "\n";
+    //!     std::abort();
+    //! }
+    //!
+ //! \param what Human readable description of the error. + //! \param where Pointer to character data where error was detected. + void parse_error_handler(const char *what, void *where); +} + +#else + +#include // For std::exception + +#define RAPIDXML_PARSE_ERROR(what, where) throw parse_error(what, where) + +namespace rapidxml +{ + + //! Parse error exception. + //! This exception is thrown by the parser when an error occurs. + //! Use what() function to get human-readable error message. + //! Use where() function to get a pointer to position within source text where error was detected. + //!

+ //! If throwing exceptions by the parser is undesirable, + //! it can be disabled by defining RAPIDXML_NO_EXCEPTIONS macro before rapidxml.hpp is included. + //! This will cause the parser to call rapidxml::parse_error_handler() function instead of throwing an exception. + //! This function must be defined by the user. + //!

+ //! This class derives from std::exception class. + class parse_error: public std::exception + { + + public: + + //! Constructs parse error + parse_error(const char *what, void *where) + : m_what(what) + , m_where(where) + { + } + + //! Gets human readable description of error. + //! \return Pointer to null terminated description of the error. + virtual const char *what() const throw() + { + return m_what; + } + + //! Gets pointer to character data where error happened. + //! Ch should be the same as char type of xml_document that produced the error. + //! \return Pointer to location within the parsed string where error occured. + template + Ch *where() const + { + return reinterpret_cast(m_where); + } + + private: + + const char *m_what; + void *m_where; + + }; +} + +#endif + +/////////////////////////////////////////////////////////////////////////// +// Pool sizes + +#ifndef RAPIDXML_STATIC_POOL_SIZE + // Size of static memory block of memory_pool. + // Define RAPIDXML_STATIC_POOL_SIZE before including rapidxml.hpp if you want to override the default value. + // No dynamic memory allocations are performed by memory_pool until static memory is exhausted. + #define RAPIDXML_STATIC_POOL_SIZE (64 * 1024) +#endif + +#ifndef RAPIDXML_DYNAMIC_POOL_SIZE + // Size of dynamic memory block of memory_pool. + // Define RAPIDXML_DYNAMIC_POOL_SIZE before including rapidxml.hpp if you want to override the default value. + // After the static block is exhausted, dynamic blocks with approximately this size are allocated by memory_pool. + #define RAPIDXML_DYNAMIC_POOL_SIZE (64 * 1024) +#endif + +#ifndef RAPIDXML_ALIGNMENT + // Memory allocation alignment. + // Define RAPIDXML_ALIGNMENT before including rapidxml.hpp if you want to override the default value, which is the size of pointer. + // All memory allocations for nodes, attributes and strings will be aligned to this value. + // This must be a power of 2 and at least 1, otherwise memory_pool will not work. + #define RAPIDXML_ALIGNMENT sizeof(void *) +#endif + +namespace rapidxml +{ + // Forward declarations + template class xml_node; + template class xml_attribute; + template class xml_document; + + //! Enumeration listing all node types produced by the parser. + //! Use xml_node::type() function to query node type. + enum node_type + { + node_document, //!< A document node. Name and value are empty. + node_element, //!< An element node. Name contains element name. Value contains text of first data node. + node_data, //!< A data node. Name is empty. Value contains data text. + node_cdata, //!< A CDATA node. Name is empty. Value contains data text. + node_comment, //!< A comment node. Name is empty. Value contains comment text. + node_declaration, //!< A declaration node. Name and value are empty. Declaration parameters (version, encoding and standalone) are in node attributes. + node_doctype, //!< A DOCTYPE node. Name is empty. Value contains DOCTYPE text. + node_pi //!< A PI node. Name contains target. Value contains instructions. + }; + + /////////////////////////////////////////////////////////////////////// + // Parsing flags + + //! Parse flag instructing the parser to not create data nodes. + //! Text of first data node will still be placed in value of parent element, unless rapidxml::parse_no_element_values flag is also specified. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_no_data_nodes = 0x1; + + //! Parse flag instructing the parser to not use text of first data node as a value of parent element. + //! Can be combined with other flags by use of | operator. + //! Note that child data nodes of element node take precendence over its value when printing. + //! That is, if element has one or more child data nodes and a value, the value will be ignored. + //! Use rapidxml::parse_no_data_nodes flag to prevent creation of data nodes if you want to manipulate data using values of elements. + //!

+ //! See xml_document::parse() function. + const int parse_no_element_values = 0x2; + + //! Parse flag instructing the parser to not place zero terminators after strings in the source text. + //! By default zero terminators are placed, modifying source text. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_no_string_terminators = 0x4; + + //! Parse flag instructing the parser to not translate entities in the source text. + //! By default entities are translated, modifying source text. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_no_entity_translation = 0x8; + + //! Parse flag instructing the parser to disable UTF-8 handling and assume plain 8 bit characters. + //! By default, UTF-8 handling is enabled. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_no_utf8 = 0x10; + + //! Parse flag instructing the parser to create XML declaration node. + //! By default, declaration node is not created. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_declaration_node = 0x20; + + //! Parse flag instructing the parser to create comments nodes. + //! By default, comment nodes are not created. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_comment_nodes = 0x40; + + //! Parse flag instructing the parser to create DOCTYPE node. + //! By default, doctype node is not created. + //! Although W3C specification allows at most one DOCTYPE node, RapidXml will silently accept documents with more than one. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_doctype_node = 0x80; + + //! Parse flag instructing the parser to create PI nodes. + //! By default, PI nodes are not created. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_pi_nodes = 0x100; + + //! Parse flag instructing the parser to validate closing tag names. + //! If not set, name inside closing tag is irrelevant to the parser. + //! By default, closing tags are not validated. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_validate_closing_tags = 0x200; + + //! Parse flag instructing the parser to trim all leading and trailing whitespace of data nodes. + //! By default, whitespace is not trimmed. + //! This flag does not cause the parser to modify source text. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_trim_whitespace = 0x400; + + //! Parse flag instructing the parser to condense all whitespace runs of data nodes to a single space character. + //! Trimming of leading and trailing whitespace of data is controlled by rapidxml::parse_trim_whitespace flag. + //! By default, whitespace is not normalized. + //! If this flag is specified, source text will be modified. + //! Can be combined with other flags by use of | operator. + //!

+ //! See xml_document::parse() function. + const int parse_normalize_whitespace = 0x800; + + // Compound flags + + //! Parse flags which represent default behaviour of the parser. + //! This is always equal to 0, so that all other flags can be simply ored together. + //! Normally there is no need to inconveniently disable flags by anding with their negated (~) values. + //! This also means that meaning of each flag is a negation of the default setting. + //! For example, if flag name is rapidxml::parse_no_utf8, it means that utf-8 is enabled by default, + //! and using the flag will disable it. + //!

+ //! See xml_document::parse() function. + const int parse_default = 0; + + //! A combination of parse flags that forbids any modifications of the source text. + //! This also results in faster parsing. However, note that the following will occur: + //!
    + //!
  • names and values of nodes will not be zero terminated, you have to use xml_base::name_size() and xml_base::value_size() functions to determine where name and value ends
    • + //!
  • entities will not be translated
    • + //!
  • whitespace will not be normalized
    • + //!
+ //! See xml_document::parse() function. + const int parse_non_destructive = parse_no_string_terminators | parse_no_entity_translation; + + //! A combination of parse flags resulting in fastest possible parsing, without sacrificing important data. + //!

+ //! See xml_document::parse() function. + const int parse_fastest = parse_non_destructive | parse_no_data_nodes; + + //! A combination of parse flags resulting in largest amount of data being extracted. + //! This usually results in slowest parsing. + //!

+ //! See xml_document::parse() function. + const int parse_full = parse_declaration_node | parse_comment_nodes | parse_doctype_node | parse_pi_nodes | parse_validate_closing_tags; + + /////////////////////////////////////////////////////////////////////// + // Internals + + //! \cond internal + namespace internal + { + + // Struct that contains lookup tables for the parser + // It must be a template to allow correct linking (because it has static data members, which are defined in a header file). + template + struct lookup_tables + { + static const unsigned char lookup_whitespace[256]; // Whitespace table + static const unsigned char lookup_node_name[256]; // Node name table + static const unsigned char lookup_text[256]; // Text table + static const unsigned char lookup_text_pure_no_ws[256]; // Text table + static const unsigned char lookup_text_pure_with_ws[256]; // Text table + static const unsigned char lookup_attribute_name[256]; // Attribute name table + static const unsigned char lookup_attribute_data_1[256]; // Attribute data table with single quote + static const unsigned char lookup_attribute_data_1_pure[256]; // Attribute data table with single quote + static const unsigned char lookup_attribute_data_2[256]; // Attribute data table with double quotes + static const unsigned char lookup_attribute_data_2_pure[256]; // Attribute data table with double quotes + static const unsigned char lookup_digits[256]; // Digits + static const unsigned char lookup_upcase[256]; // To uppercase conversion table for ASCII characters + }; + + // Find length of the string + template + inline std::size_t measure(const Ch *p) + { + const Ch *tmp = p; + while (*tmp) + ++tmp; + return tmp - p; + } + + // Compare strings for equality + template + inline bool compare(const Ch *p1, std::size_t size1, const Ch *p2, std::size_t size2, bool case_sensitive) + { + if (size1 != size2) + return false; + if (case_sensitive) + { + for (const Ch *end = p1 + size1; p1 < end; ++p1, ++p2) + if (*p1 != *p2) + return false; + } + else + { + for (const Ch *end = p1 + size1; p1 < end; ++p1, ++p2) + if (lookup_tables<0>::lookup_upcase[static_cast(*p1)] != lookup_tables<0>::lookup_upcase[static_cast(*p2)]) + return false; + } + return true; + } + } + //! \endcond + + /////////////////////////////////////////////////////////////////////// + // Memory pool + + //! This class is used by the parser to create new nodes and attributes, without overheads of dynamic memory allocation. + //! In most cases, you will not need to use this class directly. + //! However, if you need to create nodes manually or modify names/values of nodes, + //! you are encouraged to use memory_pool of relevant xml_document to allocate the memory. + //! Not only is this faster than allocating them by using new operator, + //! but also their lifetime will be tied to the lifetime of document, + //! possibly simplyfing memory management. + //!

+ //! Call allocate_node() or allocate_attribute() functions to obtain new nodes or attributes from the pool. + //! You can also call allocate_string() function to allocate strings. + //! Such strings can then be used as names or values of nodes without worrying about their lifetime. + //! Note that there is no free() function -- all allocations are freed at once when clear() function is called, + //! or when the pool is destroyed. + //!

+ //! It is also possible to create a standalone memory_pool, and use it + //! to allocate nodes, whose lifetime will not be tied to any document. + //!

+ //! Pool maintains RAPIDXML_STATIC_POOL_SIZE bytes of statically allocated memory. + //! Until static memory is exhausted, no dynamic memory allocations are done. + //! When static memory is exhausted, pool allocates additional blocks of memory of size RAPIDXML_DYNAMIC_POOL_SIZE each, + //! by using global new[] and delete[] operators. + //! This behaviour can be changed by setting custom allocation routines. + //! Use set_allocator() function to set them. + //!

+ //! Allocations for nodes, attributes and strings are aligned at RAPIDXML_ALIGNMENT bytes. + //! This value defaults to the size of pointer on target architecture. + //!

+ //! To obtain absolutely top performance from the parser, + //! it is important that all nodes are allocated from a single, contiguous block of memory. + //! Otherwise, cache misses when jumping between two (or more) disjoint blocks of memory can slow down parsing quite considerably. + //! If required, you can tweak RAPIDXML_STATIC_POOL_SIZE, RAPIDXML_DYNAMIC_POOL_SIZE and RAPIDXML_ALIGNMENT + //! to obtain best wasted memory to performance compromise. + //! To do it, define their values before rapidxml.hpp file is included. + //! \param Ch Character type of created nodes. + template + class memory_pool + { + + public: + + //! \cond internal + typedef void *(alloc_func)(std::size_t); // Type of user-defined function used to allocate memory + typedef void (free_func)(void *); // Type of user-defined function used to free memory + //! \endcond + + //! Constructs empty pool with default allocator functions. + memory_pool() + : m_alloc_func(0) + , m_free_func(0) + { + init(); + } + + //! Destroys pool and frees all the memory. + //! This causes memory occupied by nodes allocated by the pool to be freed. + //! Nodes allocated from the pool are no longer valid. + ~memory_pool() + { + clear(); + } + + //! Allocates a new node from the pool, and optionally assigns name and value to it. + //! If the allocation request cannot be accomodated, this function will throw std::bad_alloc. + //! If exceptions are disabled by defining RAPIDXML_NO_EXCEPTIONS, this function + //! will call rapidxml::parse_error_handler() function. + //! \param type Type of node to create. + //! \param name Name to assign to the node, or 0 to assign no name. + //! \param value Value to assign to the node, or 0 to assign no value. + //! \param name_size Size of name to assign, or 0 to automatically calculate size from name string. + //! \param value_size Size of value to assign, or 0 to automatically calculate size from value string. + //! \return Pointer to allocated node. This pointer will never be NULL. + xml_node *allocate_node(node_type type, + const Ch *name = 0, const Ch *value = 0, + std::size_t name_size = 0, std::size_t value_size = 0) + { + void *memory = allocate_aligned(sizeof(xml_node)); + xml_node *node = new(memory) xml_node(type); + if (name) + { + if (name_size > 0) + node->name(name, name_size); + else + node->name(name); + } + if (value) + { + if (value_size > 0) + node->value(value, value_size); + else + node->value(value); + } + return node; + } + + //! Allocates a new attribute from the pool, and optionally assigns name and value to it. + //! If the allocation request cannot be accomodated, this function will throw std::bad_alloc. + //! If exceptions are disabled by defining RAPIDXML_NO_EXCEPTIONS, this function + //! will call rapidxml::parse_error_handler() function. + //! \param name Name to assign to the attribute, or 0 to assign no name. + //! \param value Value to assign to the attribute, or 0 to assign no value. + //! \param name_size Size of name to assign, or 0 to automatically calculate size from name string. + //! \param value_size Size of value to assign, or 0 to automatically calculate size from value string. + //! \return Pointer to allocated attribute. This pointer will never be NULL. + xml_attribute *allocate_attribute(const Ch *name = 0, const Ch *value = 0, + std::size_t name_size = 0, std::size_t value_size = 0) + { + void *memory = allocate_aligned(sizeof(xml_attribute)); + xml_attribute *attribute = new(memory) xml_attribute; + if (name) + { + if (name_size > 0) + attribute->name(name, name_size); + else + attribute->name(name); + } + if (value) + { + if (value_size > 0) + attribute->value(value, value_size); + else + attribute->value(value); + } + return attribute; + } + + //! Allocates a char array of given size from the pool, and optionally copies a given string to it. + //! If the allocation request cannot be accomodated, this function will throw std::bad_alloc. + //! If exceptions are disabled by defining RAPIDXML_NO_EXCEPTIONS, this function + //! will call rapidxml::parse_error_handler() function. + //! \param source String to initialize the allocated memory with, or 0 to not initialize it. + //! \param size Number of characters to allocate, or zero to calculate it automatically from source string length; if size is 0, source string must be specified and null terminated. + //! \return Pointer to allocated char array. This pointer will never be NULL. + Ch *allocate_string(const Ch *source = 0, std::size_t size = 0) + { + assert(source || size); // Either source or size (or both) must be specified + if (size == 0) + size = internal::measure(source) + 1; + Ch *result = static_cast(allocate_aligned(size * sizeof(Ch))); + if (source) + for (std::size_t i = 0; i < size; ++i) + result[i] = source[i]; + return result; + } + + //! Clones an xml_node and its hierarchy of child nodes and attributes. + //! Nodes and attributes are allocated from this memory pool. + //! Names and values are not cloned, they are shared between the clone and the source. + //! Result node can be optionally specified as a second parameter, + //! in which case its contents will be replaced with cloned source node. + //! This is useful when you want to clone entire document. + //! \param source Node to clone. + //! \param result Node to put results in, or 0 to automatically allocate result node + //! \return Pointer to cloned node. This pointer will never be NULL. + xml_node *clone_node(const xml_node *source, xml_node *result = 0) + { + // Prepare result node + if (result) + { + result->remove_all_attributes(); + result->remove_all_nodes(); + result->type(source->type()); + } + else + result = allocate_node(source->type()); + + // Clone name and value + result->name(source->name(), source->name_size()); + result->value(source->value(), source->value_size()); + + // Clone child nodes and attributes + for (xml_node *child = source->first_node(); child; child = child->next_sibling()) + result->append_node(clone_node(child)); + for (xml_attribute *attr = source->first_attribute(); attr; attr = attr->next_attribute()) + result->append_attribute(allocate_attribute(attr->name(), attr->value(), attr->name_size(), attr->value_size())); + + return result; + } + + //! Clears the pool. + //! This causes memory occupied by nodes allocated by the pool to be freed. + //! Any nodes or strings allocated from the pool will no longer be valid. + void clear() + { + while (m_begin != m_static_memory) + { + char *previous_begin = reinterpret_cast
(align(m_begin))->previous_begin; + if (m_free_func) + m_free_func(m_begin); + else + delete[] m_begin; + m_begin = previous_begin; + } + init(); + } + + //! Sets or resets the user-defined memory allocation functions for the pool. + //! This can only be called when no memory is allocated from the pool yet, otherwise results are undefined. + //! Allocation function must not return invalid pointer on failure. It should either throw, + //! stop the program, or use longjmp() function to pass control to other place of program. + //! If it returns invalid pointer, results are undefined. + //!

+ //! User defined allocation functions must have the following forms: + //!
+ //!
void *allocate(std::size_t size); + //!
void free(void *pointer); + //!
+ //! \param af Allocation function, or 0 to restore default function + //! \param ff Free function, or 0 to restore default function + void set_allocator(alloc_func *af, free_func *ff) + { + assert(m_begin == m_static_memory && m_ptr == align(m_begin)); // Verify that no memory is allocated yet + m_alloc_func = af; + m_free_func = ff; + } + + private: + + struct header + { + char *previous_begin; + }; + + void init() + { + m_begin = m_static_memory; + m_ptr = align(m_begin); + m_end = m_static_memory + sizeof(m_static_memory); + } + + char *align(char *ptr) + { + std::size_t alignment = ((RAPIDXML_ALIGNMENT - (std::size_t(ptr) & (RAPIDXML_ALIGNMENT - 1))) & (RAPIDXML_ALIGNMENT - 1)); + return ptr + alignment; + } + + char *allocate_raw(std::size_t size) + { + // Allocate + void *memory; + if (m_alloc_func) // Allocate memory using either user-specified allocation function or global operator new[] + { + memory = m_alloc_func(size); + assert(memory); // Allocator is not allowed to return 0, on failure it must either throw, stop the program or use longjmp + } + else + { + memory = new char[size]; +#ifdef RAPIDXML_NO_EXCEPTIONS + if (!memory) // If exceptions are disabled, verify memory allocation, because new will not be able to throw bad_alloc + RAPIDXML_PARSE_ERROR("out of memory", 0); +#endif + } + return static_cast(memory); + } + + void *allocate_aligned(std::size_t size) + { + // Calculate aligned pointer + char *result = align(m_ptr); + + // If not enough memory left in current pool, allocate a new pool + if (result + size > m_end) + { + // Calculate required pool size (may be bigger than RAPIDXML_DYNAMIC_POOL_SIZE) + std::size_t pool_size = RAPIDXML_DYNAMIC_POOL_SIZE; + if (pool_size < size) + pool_size = size; + + // Allocate + std::size_t alloc_size = sizeof(header) + (2 * RAPIDXML_ALIGNMENT - 2) + pool_size; // 2 alignments required in worst case: one for header, one for actual allocation + char *raw_memory = allocate_raw(alloc_size); + + // Setup new pool in allocated memory + char *pool = align(raw_memory); + header *new_header = reinterpret_cast
(pool); + new_header->previous_begin = m_begin; + m_begin = raw_memory; + m_ptr = pool + sizeof(header); + m_end = raw_memory + alloc_size; + + // Calculate aligned pointer again using new pool + result = align(m_ptr); + } + + // Update pool and return aligned pointer + m_ptr = result + size; + return result; + } + + char *m_begin; // Start of raw memory making up current pool + char *m_ptr; // First free byte in current pool + char *m_end; // One past last available byte in current pool + char m_static_memory[RAPIDXML_STATIC_POOL_SIZE]; // Static raw memory + alloc_func *m_alloc_func; // Allocator function, or 0 if default is to be used + free_func *m_free_func; // Free function, or 0 if default is to be used + }; + + /////////////////////////////////////////////////////////////////////////// + // XML base + + //! Base class for xml_node and xml_attribute implementing common functions: + //! name(), name_size(), value(), value_size() and parent(). + //! \param Ch Character type to use + template + class xml_base + { + + public: + + /////////////////////////////////////////////////////////////////////////// + // Construction & destruction + + // Construct a base with empty name, value and parent + xml_base() + : m_name(0) + , m_value(0) + , m_parent(0) + { + } + + /////////////////////////////////////////////////////////////////////////// + // Node data access + + //! Gets name of the node. + //! Interpretation of name depends on type of node. + //! Note that name will not be zero-terminated if rapidxml::parse_no_string_terminators option was selected during parse. + //!

+ //! Use name_size() function to determine length of the name. + //! \return Name of node, or empty string if node has no name. + Ch *name() const + { + return m_name ? m_name : nullstr(); + } + + //! Gets size of node name, not including terminator character. + //! This function works correctly irrespective of whether name is or is not zero terminated. + //! \return Size of node name, in characters. + std::size_t name_size() const + { + return m_name ? m_name_size : 0; + } + + //! Gets value of node. + //! Interpretation of value depends on type of node. + //! Note that value will not be zero-terminated if rapidxml::parse_no_string_terminators option was selected during parse. + //!

+ //! Use value_size() function to determine length of the value. + //! \return Value of node, or empty string if node has no value. + Ch *value() const + { + return m_value ? m_value : nullstr(); + } + + //! Gets size of node value, not including terminator character. + //! This function works correctly irrespective of whether value is or is not zero terminated. + //! \return Size of node value, in characters. + std::size_t value_size() const + { + return m_value ? m_value_size : 0; + } + + /////////////////////////////////////////////////////////////////////////// + // Node modification + + //! Sets name of node to a non zero-terminated string. + //! See \ref ownership_of_strings. + //!

+ //! Note that node does not own its name or value, it only stores a pointer to it. + //! It will not delete or otherwise free the pointer on destruction. + //! It is reponsibility of the user to properly manage lifetime of the string. + //! The easiest way to achieve it is to use memory_pool of the document to allocate the string - + //! on destruction of the document the string will be automatically freed. + //!

+ //! Size of name must be specified separately, because name does not have to be zero terminated. + //! Use name(const Ch *) function to have the length automatically calculated (string must be zero terminated). + //! \param name Name of node to set. Does not have to be zero terminated. + //! \param size Size of name, in characters. This does not include zero terminator, if one is present. + void name(const Ch *name, std::size_t size) + { + m_name = const_cast(name); + m_name_size = size; + } + + //! Sets name of node to a zero-terminated string. + //! See also \ref ownership_of_strings and xml_node::name(const Ch *, std::size_t). + //! \param name Name of node to set. Must be zero terminated. + void name(const Ch *name) + { + this->name(name, internal::measure(name)); + } + + //! Sets value of node to a non zero-terminated string. + //! See \ref ownership_of_strings. + //!

+ //! Note that node does not own its name or value, it only stores a pointer to it. + //! It will not delete or otherwise free the pointer on destruction. + //! It is reponsibility of the user to properly manage lifetime of the string. + //! The easiest way to achieve it is to use memory_pool of the document to allocate the string - + //! on destruction of the document the string will be automatically freed. + //!

+ //! Size of value must be specified separately, because it does not have to be zero terminated. + //! Use value(const Ch *) function to have the length automatically calculated (string must be zero terminated). + //!

+ //! If an element has a child node of type node_data, it will take precedence over element value when printing. + //! If you want to manipulate data of elements using values, use parser flag rapidxml::parse_no_data_nodes to prevent creation of data nodes by the parser. + //! \param value value of node to set. Does not have to be zero terminated. + //! \param size Size of value, in characters. This does not include zero terminator, if one is present. + void value(const Ch *value, std::size_t size) + { + m_value = const_cast(value); + m_value_size = size; + } + + //! Sets value of node to a zero-terminated string. + //! See also \ref ownership_of_strings and xml_node::value(const Ch *, std::size_t). + //! \param value Vame of node to set. Must be zero terminated. + void value(const Ch *value) + { + this->value(value, internal::measure(value)); + } + + /////////////////////////////////////////////////////////////////////////// + // Related nodes access + + //! Gets node parent. + //! \return Pointer to parent node, or 0 if there is no parent. + xml_node *parent() const + { + return m_parent; + } + + protected: + + // Return empty string + static Ch *nullstr() + { + static Ch zero = Ch('\0'); + return &zero; + } + + Ch *m_name; // Name of node, or 0 if no name + Ch *m_value; // Value of node, or 0 if no value + std::size_t m_name_size; // Length of node name, or undefined of no name + std::size_t m_value_size; // Length of node value, or undefined if no value + xml_node *m_parent; // Pointer to parent node, or 0 if none + + }; + + //! Class representing attribute node of XML document. + //! Each attribute has name and value strings, which are available through name() and value() functions (inherited from xml_base). + //! Note that after parse, both name and value of attribute will point to interior of source text used for parsing. + //! Thus, this text must persist in memory for the lifetime of attribute. + //! \param Ch Character type to use. + template + class xml_attribute: public xml_base + { + + friend class xml_node; + + public: + + /////////////////////////////////////////////////////////////////////////// + // Construction & destruction + + //! Constructs an empty attribute with the specified type. + //! Consider using memory_pool of appropriate xml_document if allocating attributes manually. + xml_attribute() + { + } + + /////////////////////////////////////////////////////////////////////////// + // Related nodes access + + //! Gets document of which attribute is a child. + //! \return Pointer to document that contains this attribute, or 0 if there is no parent document. + xml_document *document() const + { + if (xml_node *node = this->parent()) + { + while (node->parent()) + node = node->parent(); + return node->type() == node_document ? static_cast *>(node) : 0; + } + else + return 0; + } + + //! Gets previous attribute, optionally matching attribute name. + //! \param name Name of attribute to find, or 0 to return previous attribute regardless of its name; this string doesn't have to be zero-terminated if name_size is non-zero + //! \param name_size Size of name, in characters, or 0 to have size calculated automatically from string + //! \param case_sensitive Should name comparison be case-sensitive; non case-sensitive comparison works properly only for ASCII characters + //! \return Pointer to found attribute, or 0 if not found. + xml_attribute *previous_attribute(const Ch *name = 0, std::size_t name_size = 0, bool case_sensitive = true) const + { + if (name) + { + if (name_size == 0) + name_size = internal::measure(name); + for (xml_attribute *attribute = m_prev_attribute; attribute; attribute = attribute->m_prev_attribute) + if (internal::compare(attribute->name(), attribute->name_size(), name, name_size, case_sensitive)) + return attribute; + return 0; + } + else + return this->m_parent ? m_prev_attribute : 0; + } + + //! Gets next attribute, optionally matching attribute name. + //! \param name Name of attribute to find, or 0 to return next attribute regardless of its name; this string doesn't have to be zero-terminated if name_size is non-zero + //! \param name_size Size of name, in characters, or 0 to have size calculated automatically from string + //! \param case_sensitive Should name comparison be case-sensitive; non case-sensitive comparison works properly only for ASCII characters + //! \return Pointer to found attribute, or 0 if not found. + xml_attribute *next_attribute(const Ch *name = 0, std::size_t name_size = 0, bool case_sensitive = true) const + { + if (name) + { + if (name_size == 0) + name_size = internal::measure(name); + for (xml_attribute *attribute = m_next_attribute; attribute; attribute = attribute->m_next_attribute) + if (internal::compare(attribute->name(), attribute->name_size(), name, name_size, case_sensitive)) + return attribute; + return 0; + } + else + return this->m_parent ? m_next_attribute : 0; + } + + private: + + xml_attribute *m_prev_attribute; // Pointer to previous sibling of attribute, or 0 if none; only valid if parent is non-zero + xml_attribute *m_next_attribute; // Pointer to next sibling of attribute, or 0 if none; only valid if parent is non-zero + + }; + + /////////////////////////////////////////////////////////////////////////// + // XML node + + //! Class representing a node of XML document. + //! Each node may have associated name and value strings, which are available through name() and value() functions. + //! Interpretation of name and value depends on type of the node. + //! Type of node can be determined by using type() function. + //!

+ //! Note that after parse, both name and value of node, if any, will point interior of source text used for parsing. + //! Thus, this text must persist in the memory for the lifetime of node. + //! \param Ch Character type to use. + template + class xml_node: public xml_base + { + + public: + + /////////////////////////////////////////////////////////////////////////// + // Construction & destruction + + //! Constructs an empty node with the specified type. + //! Consider using memory_pool of appropriate document to allocate nodes manually. + //! \param type Type of node to construct. + xml_node(node_type type) + : m_type(type) + , m_first_node(0) + , m_first_attribute(0) + { + } + + /////////////////////////////////////////////////////////////////////////// + // Node data access + + //! Gets type of node. + //! \return Type of node. + node_type type() const + { + return m_type; + } + + /////////////////////////////////////////////////////////////////////////// + // Related nodes access + + //! Gets document of which node is a child. + //! \return Pointer to document that contains this node, or 0 if there is no parent document. + xml_document *document() const + { + xml_node *node = const_cast *>(this); + while (node->parent()) + node = node->parent(); + return node->type() == node_document ? static_cast *>(node) : 0; + } + + //! Gets first child node, optionally matching node name. + //! \param name Name of child to find, or 0 to return first child regardless of its name; this string doesn't have to be zero-terminated if name_size is non-zero + //! \param name_size Size of name, in characters, or 0 to have size calculated automatically from string + //! \param case_sensitive Should name comparison be case-sensitive; non case-sensitive comparison works properly only for ASCII characters + //! \return Pointer to found child, or 0 if not found. + xml_node *first_node(const Ch *name = 0, std::size_t name_size = 0, bool case_sensitive = true) const + { + if (name) + { + if (name_size == 0) + name_size = internal::measure(name); + for (xml_node *child = m_first_node; child; child = child->next_sibling()) + if (internal::compare(child->name(), child->name_size(), name, name_size, case_sensitive)) + return child; + return 0; + } + else + return m_first_node; + } + + //! Gets last child node, optionally matching node name. + //! Behaviour is undefined if node has no children. + //! Use first_node() to test if node has children. + //! \param name Name of child to find, or 0 to return last child regardless of its name; this string doesn't have to be zero-terminated if name_size is non-zero + //! \param name_size Size of name, in characters, or 0 to have size calculated automatically from string + //! \param case_sensitive Should name comparison be case-sensitive; non case-sensitive comparison works properly only for ASCII characters + //! \return Pointer to found child, or 0 if not found. + xml_node *last_node(const Ch *name = 0, std::size_t name_size = 0, bool case_sensitive = true) const + { + assert(m_first_node); // Cannot query for last child if node has no children + if (name) + { + if (name_size == 0) + name_size = internal::measure(name); + for (xml_node *child = m_last_node; child; child = child->previous_sibling()) + if (internal::compare(child->name(), child->name_size(), name, name_size, case_sensitive)) + return child; + return 0; + } + else + return m_last_node; + } + + //! Gets previous sibling node, optionally matching node name. + //! Behaviour is undefined if node has no parent. + //! Use parent() to test if node has a parent. + //! \param name Name of sibling to find, or 0 to return previous sibling regardless of its name; this string doesn't have to be zero-terminated if name_size is non-zero + //! \param name_size Size of name, in characters, or 0 to have size calculated automatically from string + //! \param case_sensitive Should name comparison be case-sensitive; non case-sensitive comparison works properly only for ASCII characters + //! \return Pointer to found sibling, or 0 if not found. + xml_node *previous_sibling(const Ch *name = 0, std::size_t name_size = 0, bool case_sensitive = true) const + { + assert(this->m_parent); // Cannot query for siblings if node has no parent + if (name) + { + if (name_size == 0) + name_size = internal::measure(name); + for (xml_node *sibling = m_prev_sibling; sibling; sibling = sibling->m_prev_sibling) + if (internal::compare(sibling->name(), sibling->name_size(), name, name_size, case_sensitive)) + return sibling; + return 0; + } + else + return m_prev_sibling; + } + + //! Gets next sibling node, optionally matching node name. + //! Behaviour is undefined if node has no parent. + //! Use parent() to test if node has a parent. + //! \param name Name of sibling to find, or 0 to return next sibling regardless of its name; this string doesn't have to be zero-terminated if name_size is non-zero + //! \param name_size Size of name, in characters, or 0 to have size calculated automatically from string + //! \param case_sensitive Should name comparison be case-sensitive; non case-sensitive comparison works properly only for ASCII characters + //! \return Pointer to found sibling, or 0 if not found. + xml_node *next_sibling(const Ch *name = 0, std::size_t name_size = 0, bool case_sensitive = true) const + { + assert(this->m_parent); // Cannot query for siblings if node has no parent + if (name) + { + if (name_size == 0) + name_size = internal::measure(name); + for (xml_node *sibling = m_next_sibling; sibling; sibling = sibling->m_next_sibling) + if (internal::compare(sibling->name(), sibling->name_size(), name, name_size, case_sensitive)) + return sibling; + return 0; + } + else + return m_next_sibling; + } + + //! Gets first attribute of node, optionally matching attribute name. + //! \param name Name of attribute to find, or 0 to return first attribute regardless of its name; this string doesn't have to be zero-terminated if name_size is non-zero + //! \param name_size Size of name, in characters, or 0 to have size calculated automatically from string + //! \param case_sensitive Should name comparison be case-sensitive; non case-sensitive comparison works properly only for ASCII characters + //! \return Pointer to found attribute, or 0 if not found. + xml_attribute *first_attribute(const Ch *name = 0, std::size_t name_size = 0, bool case_sensitive = true) const + { + if (name) + { + if (name_size == 0) + name_size = internal::measure(name); + for (xml_attribute *attribute = m_first_attribute; attribute; attribute = attribute->m_next_attribute) + if (internal::compare(attribute->name(), attribute->name_size(), name, name_size, case_sensitive)) + return attribute; + return 0; + } + else + return m_first_attribute; + } + + //! Gets last attribute of node, optionally matching attribute name. + //! \param name Name of attribute to find, or 0 to return last attribute regardless of its name; this string doesn't have to be zero-terminated if name_size is non-zero + //! \param name_size Size of name, in characters, or 0 to have size calculated automatically from string + //! \param case_sensitive Should name comparison be case-sensitive; non case-sensitive comparison works properly only for ASCII characters + //! \return Pointer to found attribute, or 0 if not found. + xml_attribute *last_attribute(const Ch *name = 0, std::size_t name_size = 0, bool case_sensitive = true) const + { + if (name) + { + if (name_size == 0) + name_size = internal::measure(name); + for (xml_attribute *attribute = m_last_attribute; attribute; attribute = attribute->m_prev_attribute) + if (internal::compare(attribute->name(), attribute->name_size(), name, name_size, case_sensitive)) + return attribute; + return 0; + } + else + return m_first_attribute ? m_last_attribute : 0; + } + + /////////////////////////////////////////////////////////////////////////// + // Node modification + + //! Sets type of node. + //! \param type Type of node to set. + void type(node_type type) + { + m_type = type; + } + + /////////////////////////////////////////////////////////////////////////// + // Node manipulation + + //! Prepends a new child node. + //! The prepended child becomes the first child, and all existing children are moved one position back. + //! \param child Node to prepend. + void prepend_node(xml_node *child) + { + assert(child && !child->parent() && child->type() != node_document); + if (first_node()) + { + child->m_next_sibling = m_first_node; + m_first_node->m_prev_sibling = child; + } + else + { + child->m_next_sibling = 0; + m_last_node = child; + } + m_first_node = child; + child->m_parent = this; + child->m_prev_sibling = 0; + } + + //! Appends a new child node. + //! The appended child becomes the last child. + //! \param child Node to append. + void append_node(xml_node *child) + { + assert(child && !child->parent() && child->type() != node_document); + if (first_node()) + { + child->m_prev_sibling = m_last_node; + m_last_node->m_next_sibling = child; + } + else + { + child->m_prev_sibling = 0; + m_first_node = child; + } + m_last_node = child; + child->m_parent = this; + child->m_next_sibling = 0; + } + + //! Inserts a new child node at specified place inside the node. + //! All children after and including the specified node are moved one position back. + //! \param where Place where to insert the child, or 0 to insert at the back. + //! \param child Node to insert. + void insert_node(xml_node *where, xml_node *child) + { + assert(!where || where->parent() == this); + assert(child && !child->parent() && child->type() != node_document); + if (where == m_first_node) + prepend_node(child); + else if (where == 0) + append_node(child); + else + { + child->m_prev_sibling = where->m_prev_sibling; + child->m_next_sibling = where; + where->m_prev_sibling->m_next_sibling = child; + where->m_prev_sibling = child; + child->m_parent = this; + } + } + + //! Removes first child node. + //! If node has no children, behaviour is undefined. + //! Use first_node() to test if node has children. + void remove_first_node() + { + assert(first_node()); + xml_node *child = m_first_node; + m_first_node = child->m_next_sibling; + if (child->m_next_sibling) + child->m_next_sibling->m_prev_sibling = 0; + else + m_last_node = 0; + child->m_parent = 0; + } + + //! Removes last child of the node. + //! If node has no children, behaviour is undefined. + //! Use first_node() to test if node has children. + void remove_last_node() + { + assert(first_node()); + xml_node *child = m_last_node; + if (child->m_prev_sibling) + { + m_last_node = child->m_prev_sibling; + child->m_prev_sibling->m_next_sibling = 0; + } + else + m_first_node = 0; + child->m_parent = 0; + } + + //! Removes specified child from the node + // \param where Pointer to child to be removed. + void remove_node(xml_node *where) + { + assert(where && where->parent() == this); + assert(first_node()); + if (where == m_first_node) + remove_first_node(); + else if (where == m_last_node) + remove_last_node(); + else + { + where->m_prev_sibling->m_next_sibling = where->m_next_sibling; + where->m_next_sibling->m_prev_sibling = where->m_prev_sibling; + where->m_parent = 0; + } + } + + //! Removes all child nodes (but not attributes). + void remove_all_nodes() + { + for (xml_node *node = first_node(); node; node = node->m_next_sibling) + node->m_parent = 0; + m_first_node = 0; + } + + //! Prepends a new attribute to the node. + //! \param attribute Attribute to prepend. + void prepend_attribute(xml_attribute *attribute) + { + assert(attribute && !attribute->parent()); + if (first_attribute()) + { + attribute->m_next_attribute = m_first_attribute; + m_first_attribute->m_prev_attribute = attribute; + } + else + { + attribute->m_next_attribute = 0; + m_last_attribute = attribute; + } + m_first_attribute = attribute; + attribute->m_parent = this; + attribute->m_prev_attribute = 0; + } + + //! Appends a new attribute to the node. + //! \param attribute Attribute to append. + void append_attribute(xml_attribute *attribute) + { + assert(attribute && !attribute->parent()); + if (first_attribute()) + { + attribute->m_prev_attribute = m_last_attribute; + m_last_attribute->m_next_attribute = attribute; + } + else + { + attribute->m_prev_attribute = 0; + m_first_attribute = attribute; + } + m_last_attribute = attribute; + attribute->m_parent = this; + attribute->m_next_attribute = 0; + } + + //! Inserts a new attribute at specified place inside the node. + //! All attributes after and including the specified attribute are moved one position back. + //! \param where Place where to insert the attribute, or 0 to insert at the back. + //! \param attribute Attribute to insert. + void insert_attribute(xml_attribute *where, xml_attribute *attribute) + { + assert(!where || where->parent() == this); + assert(attribute && !attribute->parent()); + if (where == m_first_attribute) + prepend_attribute(attribute); + else if (where == 0) + append_attribute(attribute); + else + { + attribute->m_prev_attribute = where->m_prev_attribute; + attribute->m_next_attribute = where; + where->m_prev_attribute->m_next_attribute = attribute; + where->m_prev_attribute = attribute; + attribute->m_parent = this; + } + } + + //! Removes first attribute of the node. + //! If node has no attributes, behaviour is undefined. + //! Use first_attribute() to test if node has attributes. + void remove_first_attribute() + { + assert(first_attribute()); + xml_attribute *attribute = m_first_attribute; + if (attribute->m_next_attribute) + { + attribute->m_next_attribute->m_prev_attribute = 0; + } + else + m_last_attribute = 0; + attribute->m_parent = 0; + m_first_attribute = attribute->m_next_attribute; + } + + //! Removes last attribute of the node. + //! If node has no attributes, behaviour is undefined. + //! Use first_attribute() to test if node has attributes. + void remove_last_attribute() + { + assert(first_attribute()); + xml_attribute *attribute = m_last_attribute; + if (attribute->m_prev_attribute) + { + attribute->m_prev_attribute->m_next_attribute = 0; + m_last_attribute = attribute->m_prev_attribute; + } + else + m_first_attribute = 0; + attribute->m_parent = 0; + } + + //! Removes specified attribute from node. + //! \param where Pointer to attribute to be removed. + void remove_attribute(xml_attribute *where) + { + assert(first_attribute() && where->parent() == this); + if (where == m_first_attribute) + remove_first_attribute(); + else if (where == m_last_attribute) + remove_last_attribute(); + else + { + where->m_prev_attribute->m_next_attribute = where->m_next_attribute; + where->m_next_attribute->m_prev_attribute = where->m_prev_attribute; + where->m_parent = 0; + } + } + + //! Removes all attributes of node. + void remove_all_attributes() + { + for (xml_attribute *attribute = first_attribute(); attribute; attribute = attribute->m_next_attribute) + attribute->m_parent = 0; + m_first_attribute = 0; + } + + private: + + /////////////////////////////////////////////////////////////////////////// + // Restrictions + + // No copying + xml_node(const xml_node &); + void operator =(const xml_node &); + + /////////////////////////////////////////////////////////////////////////// + // Data members + + // Note that some of the pointers below have UNDEFINED values if certain other pointers are 0. + // This is required for maximum performance, as it allows the parser to omit initialization of + // unneded/redundant values. + // + // The rules are as follows: + // 1. first_node and first_attribute contain valid pointers, or 0 if node has no children/attributes respectively + // 2. last_node and last_attribute are valid only if node has at least one child/attribute respectively, otherwise they contain garbage + // 3. prev_sibling and next_sibling are valid only if node has a parent, otherwise they contain garbage + + node_type m_type; // Type of node; always valid + xml_node *m_first_node; // Pointer to first child node, or 0 if none; always valid + xml_node *m_last_node; // Pointer to last child node, or 0 if none; this value is only valid if m_first_node is non-zero + xml_attribute *m_first_attribute; // Pointer to first attribute of node, or 0 if none; always valid + xml_attribute *m_last_attribute; // Pointer to last attribute of node, or 0 if none; this value is only valid if m_first_attribute is non-zero + xml_node *m_prev_sibling; // Pointer to previous sibling of node, or 0 if none; this value is only valid if m_parent is non-zero + xml_node *m_next_sibling; // Pointer to next sibling of node, or 0 if none; this value is only valid if m_parent is non-zero + + }; + + /////////////////////////////////////////////////////////////////////////// + // XML document + + //! This class represents root of the DOM hierarchy. + //! It is also an xml_node and a memory_pool through public inheritance. + //! Use parse() function to build a DOM tree from a zero-terminated XML text string. + //! parse() function allocates memory for nodes and attributes by using functions of xml_document, + //! which are inherited from memory_pool. + //! To access root node of the document, use the document itself, as if it was an xml_node. + //! \param Ch Character type to use. + template + class xml_document: public xml_node, public memory_pool + { + + public: + + //! Constructs empty XML document + xml_document() + : xml_node(node_document) + { + } + + //! Parses zero-terminated XML string according to given flags. + //! Passed string will be modified by the parser, unless rapidxml::parse_non_destructive flag is used. + //! The string must persist for the lifetime of the document. + //! In case of error, rapidxml::parse_error exception will be thrown. + //!

+ //! If you want to parse contents of a file, you must first load the file into the memory, and pass pointer to its beginning. + //! Make sure that data is zero-terminated. + //!

+ //! Document can be parsed into multiple times. + //! Each new call to parse removes previous nodes and attributes (if any), but does not clear memory pool. + //! \param text XML data to parse; pointer is non-const to denote fact that this data may be modified by the parser. + template + void parse(Ch *text) + { + assert(text); + + // Remove current contents + this->remove_all_nodes(); + this->remove_all_attributes(); + + // Parse BOM, if any + parse_bom(text); + + // Parse children + while (1) + { + // Skip whitespace before node + skip(text); + if (*text == 0) + break; + + // Parse and append new child + if (*text == Ch('<')) + { + ++text; // Skip '<' + if (xml_node *node = parse_node(text)) + this->append_node(node); + } + else + RAPIDXML_PARSE_ERROR("expected <", text); + } + + } + + //! Clears the document by deleting all nodes and clearing the memory pool. + //! All nodes owned by document pool are destroyed. + void clear() + { + this->remove_all_nodes(); + this->remove_all_attributes(); + memory_pool::clear(); + } + + private: + + /////////////////////////////////////////////////////////////////////// + // Internal character utility functions + + // Detect whitespace character + struct whitespace_pred + { + static unsigned char test(Ch ch) + { + return internal::lookup_tables<0>::lookup_whitespace[static_cast(ch)]; + } + }; + + // Detect node name character + struct node_name_pred + { + static unsigned char test(Ch ch) + { + return internal::lookup_tables<0>::lookup_node_name[static_cast(ch)]; + } + }; + + // Detect attribute name character + struct attribute_name_pred + { + static unsigned char test(Ch ch) + { + return internal::lookup_tables<0>::lookup_attribute_name[static_cast(ch)]; + } + }; + + // Detect text character (PCDATA) + struct text_pred + { + static unsigned char test(Ch ch) + { + return internal::lookup_tables<0>::lookup_text[static_cast(ch)]; + } + }; + + // Detect text character (PCDATA) that does not require processing + struct text_pure_no_ws_pred + { + static unsigned char test(Ch ch) + { + return internal::lookup_tables<0>::lookup_text_pure_no_ws[static_cast(ch)]; + } + }; + + // Detect text character (PCDATA) that does not require processing + struct text_pure_with_ws_pred + { + static unsigned char test(Ch ch) + { + return internal::lookup_tables<0>::lookup_text_pure_with_ws[static_cast(ch)]; + } + }; + + // Detect attribute value character + template + struct attribute_value_pred + { + static unsigned char test(Ch ch) + { + if (Quote == Ch('\'')) + return internal::lookup_tables<0>::lookup_attribute_data_1[static_cast(ch)]; + if (Quote == Ch('\"')) + return internal::lookup_tables<0>::lookup_attribute_data_2[static_cast(ch)]; + return 0; // Should never be executed, to avoid warnings on Comeau + } + }; + + // Detect attribute value character + template + struct attribute_value_pure_pred + { + static unsigned char test(Ch ch) + { + if (Quote == Ch('\'')) + return internal::lookup_tables<0>::lookup_attribute_data_1_pure[static_cast(ch)]; + if (Quote == Ch('\"')) + return internal::lookup_tables<0>::lookup_attribute_data_2_pure[static_cast(ch)]; + return 0; // Should never be executed, to avoid warnings on Comeau + } + }; + + // Insert coded character, using UTF8 or 8-bit ASCII + template + static void insert_coded_character(Ch *&text, unsigned long code) + { + if (Flags & parse_no_utf8) + { + // Insert 8-bit ASCII character + // Todo: possibly verify that code is less than 256 and use replacement char otherwise? + text[0] = static_cast(code); + text += 1; + } + else + { + // Insert UTF8 sequence + if (code < 0x80) // 1 byte sequence + { + text[0] = static_cast(code); + text += 1; + } + else if (code < 0x800) // 2 byte sequence + { + text[1] = static_cast((code | 0x80) & 0xBF); code >>= 6; + text[0] = static_cast(code | 0xC0); + text += 2; + } + else if (code < 0x10000) // 3 byte sequence + { + text[2] = static_cast((code | 0x80) & 0xBF); code >>= 6; + text[1] = static_cast((code | 0x80) & 0xBF); code >>= 6; + text[0] = static_cast(code | 0xE0); + text += 3; + } + else if (code < 0x110000) // 4 byte sequence + { + text[3] = static_cast((code | 0x80) & 0xBF); code >>= 6; + text[2] = static_cast((code | 0x80) & 0xBF); code >>= 6; + text[1] = static_cast((code | 0x80) & 0xBF); code >>= 6; + text[0] = static_cast(code | 0xF0); + text += 4; + } + else // Invalid, only codes up to 0x10FFFF are allowed in Unicode + { + RAPIDXML_PARSE_ERROR("invalid numeric character entity", text); + } + } + } + + // Skip characters until predicate evaluates to true + template + static void skip(Ch *&text) + { + Ch *tmp = text; + while (StopPred::test(*tmp)) + ++tmp; + text = tmp; + } + + // Skip characters until predicate evaluates to true while doing the following: + // - replacing XML character entity references with proper characters (' & " < > &#...;) + // - condensing whitespace sequences to single space character + template + static Ch *skip_and_expand_character_refs(Ch *&text) + { + // If entity translation, whitespace condense and whitespace trimming is disabled, use plain skip + if (Flags & parse_no_entity_translation && + !(Flags & parse_normalize_whitespace) && + !(Flags & parse_trim_whitespace)) + { + skip(text); + return text; + } + + // Use simple skip until first modification is detected + skip(text); + + // Use translation skip + Ch *src = text; + Ch *dest = src; + while (StopPred::test(*src)) + { + // If entity translation is enabled + if (!(Flags & parse_no_entity_translation)) + { + // Test if replacement is needed + if (src[0] == Ch('&')) + { + switch (src[1]) + { + + // & ' + case Ch('a'): + if (src[2] == Ch('m') && src[3] == Ch('p') && src[4] == Ch(';')) + { + *dest = Ch('&'); + ++dest; + src += 5; + continue; + } + if (src[2] == Ch('p') && src[3] == Ch('o') && src[4] == Ch('s') && src[5] == Ch(';')) + { + *dest = Ch('\''); + ++dest; + src += 6; + continue; + } + break; + + // " + case Ch('q'): + if (src[2] == Ch('u') && src[3] == Ch('o') && src[4] == Ch('t') && src[5] == Ch(';')) + { + *dest = Ch('"'); + ++dest; + src += 6; + continue; + } + break; + + // > + case Ch('g'): + if (src[2] == Ch('t') && src[3] == Ch(';')) + { + *dest = Ch('>'); + ++dest; + src += 4; + continue; + } + break; + + // < + case Ch('l'): + if (src[2] == Ch('t') && src[3] == Ch(';')) + { + *dest = Ch('<'); + ++dest; + src += 4; + continue; + } + break; + + // &#...; - assumes ASCII + case Ch('#'): + if (src[2] == Ch('x')) + { + unsigned long code = 0; + src += 3; // Skip &#x + while (1) + { + unsigned char digit = internal::lookup_tables<0>::lookup_digits[static_cast(*src)]; + if (digit == 0xFF) + break; + code = code * 16 + digit; + ++src; + } + insert_coded_character(dest, code); // Put character in output + } + else + { + unsigned long code = 0; + src += 2; // Skip &# + while (1) + { + unsigned char digit = internal::lookup_tables<0>::lookup_digits[static_cast(*src)]; + if (digit == 0xFF) + break; + code = code * 10 + digit; + ++src; + } + insert_coded_character(dest, code); // Put character in output + } + if (*src == Ch(';')) + ++src; + else + RAPIDXML_PARSE_ERROR("expected ;", src); + continue; + + // Something else + default: + // Ignore, just copy '&' verbatim + break; + + } + } + } + + // If whitespace condensing is enabled + if (Flags & parse_normalize_whitespace) + { + // Test if condensing is needed + if (whitespace_pred::test(*src)) + { + *dest = Ch(' '); ++dest; // Put single space in dest + ++src; // Skip first whitespace char + // Skip remaining whitespace chars + while (whitespace_pred::test(*src)) + ++src; + continue; + } + } + + // No replacement, only copy character + *dest++ = *src++; + + } + + // Return new end + text = src; + return dest; + + } + + /////////////////////////////////////////////////////////////////////// + // Internal parsing functions + + // Parse BOM, if any + template + void parse_bom(Ch *&text) + { + // UTF-8? + if (static_cast(text[0]) == 0xEF && + static_cast(text[1]) == 0xBB && + static_cast(text[2]) == 0xBF) + { + text += 3; // Skup utf-8 bom + } + } + + // Parse XML declaration ( + xml_node *parse_xml_declaration(Ch *&text) + { + // If parsing of declaration is disabled + if (!(Flags & parse_declaration_node)) + { + // Skip until end of declaration + while (text[0] != Ch('?') || text[1] != Ch('>')) + { + if (!text[0]) + RAPIDXML_PARSE_ERROR("unexpected end of data", text); + ++text; + } + text += 2; // Skip '?>' + return 0; + } + + // Create declaration + xml_node *declaration = this->allocate_node(node_declaration); + + // Skip whitespace before attributes or ?> + skip(text); + + // Parse declaration attributes + parse_node_attributes(text, declaration); + + // Skip ?> + if (text[0] != Ch('?') || text[1] != Ch('>')) + RAPIDXML_PARSE_ERROR("expected ?>", text); + text += 2; + + return declaration; + } + + // Parse XML comment (' + return 0; // Do not produce comment node + } + + // Remember value start + Ch *value = text; + + // Skip until end of comment + while (text[0] != Ch('-') || text[1] != Ch('-') || text[2] != Ch('>')) + { + if (!text[0]) + RAPIDXML_PARSE_ERROR("unexpected end of data", text); + ++text; + } + + // Create comment node + xml_node *comment = this->allocate_node(node_comment); + comment->value(value, text - value); + + // Place zero terminator after comment value + if (!(Flags & parse_no_string_terminators)) + *text = Ch('\0'); + + text += 3; // Skip '-->' + return comment; + } + + // Parse DOCTYPE + template + xml_node *parse_doctype(Ch *&text) + { + // Remember value start + Ch *value = text; + + // Skip to > + while (*text != Ch('>')) + { + // Determine character type + switch (*text) + { + + // If '[' encountered, scan for matching ending ']' using naive algorithm with depth + // This works for all W3C test files except for 2 most wicked + case Ch('['): + { + ++text; // Skip '[' + int depth = 1; + while (depth > 0) + { + switch (*text) + { + case Ch('['): ++depth; break; + case Ch(']'): --depth; break; + case 0: RAPIDXML_PARSE_ERROR("unexpected end of data", text); + } + ++text; + } + break; + } + + // Error on end of text + case Ch('\0'): + RAPIDXML_PARSE_ERROR("unexpected end of data", text); + + // Other character, skip it + default: + ++text; + + } + } + + // If DOCTYPE nodes enabled + if (Flags & parse_doctype_node) + { + // Create a new doctype node + xml_node *doctype = this->allocate_node(node_doctype); + doctype->value(value, text - value); + + // Place zero terminator after value + if (!(Flags & parse_no_string_terminators)) + *text = Ch('\0'); + + text += 1; // skip '>' + return doctype; + } + else + { + text += 1; // skip '>' + return 0; + } + + } + + // Parse PI + template + xml_node *parse_pi(Ch *&text) + { + // If creation of PI nodes is enabled + if (Flags & parse_pi_nodes) + { + // Create pi node + xml_node *pi = this->allocate_node(node_pi); + + // Extract PI target name + Ch *name = text; + skip(text); + if (text == name) + RAPIDXML_PARSE_ERROR("expected PI target", text); + pi->name(name, text - name); + + // Skip whitespace between pi target and pi + skip(text); + + // Remember start of pi + Ch *value = text; + + // Skip to '?>' + while (text[0] != Ch('?') || text[1] != Ch('>')) + { + if (*text == Ch('\0')) + RAPIDXML_PARSE_ERROR("unexpected end of data", text); + ++text; + } + + // Set pi value (verbatim, no entity expansion or whitespace normalization) + pi->value(value, text - value); + + // Place zero terminator after name and value + if (!(Flags & parse_no_string_terminators)) + { + pi->name()[pi->name_size()] = Ch('\0'); + pi->value()[pi->value_size()] = Ch('\0'); + } + + text += 2; // Skip '?>' + return pi; + } + else + { + // Skip to '?>' + while (text[0] != Ch('?') || text[1] != Ch('>')) + { + if (*text == Ch('\0')) + RAPIDXML_PARSE_ERROR("unexpected end of data", text); + ++text; + } + text += 2; // Skip '?>' + return 0; + } + } + + // Parse and append data + // Return character that ends data. + // This is necessary because this character might have been overwritten by a terminating 0 + template + Ch parse_and_append_data(xml_node *node, Ch *&text, Ch *contents_start) + { + // Backup to contents start if whitespace trimming is disabled + if (!(Flags & parse_trim_whitespace)) + text = contents_start; + + // Skip until end of data + Ch *value = text, *end; + if (Flags & parse_normalize_whitespace) + end = skip_and_expand_character_refs(text); + else + end = skip_and_expand_character_refs(text); + + // Trim trailing whitespace if flag is set; leading was already trimmed by whitespace skip after > + if (Flags & parse_trim_whitespace) + { + if (Flags & parse_normalize_whitespace) + { + // Whitespace is already condensed to single space characters by skipping function, so just trim 1 char off the end + if (*(end - 1) == Ch(' ')) + --end; + } + else + { + // Backup until non-whitespace character is found + while (whitespace_pred::test(*(end - 1))) + --end; + } + } + + // If characters are still left between end and value (this test is only necessary if normalization is enabled) + // Create new data node + if (!(Flags & parse_no_data_nodes)) + { + xml_node *data = this->allocate_node(node_data); + data->value(value, end - value); + node->append_node(data); + } + + // Add data to parent node if no data exists yet + if (!(Flags & parse_no_element_values)) + if (*node->value() == Ch('\0')) + node->value(value, end - value); + + // Place zero terminator after value + if (!(Flags & parse_no_string_terminators)) + { + Ch ch = *text; + *end = Ch('\0'); + return ch; // Return character that ends data; this is required because zero terminator overwritten it + } + + // Return character that ends data + return *text; + } + + // Parse CDATA + template + xml_node *parse_cdata(Ch *&text) + { + // If CDATA is disabled + if (Flags & parse_no_data_nodes) + { + // Skip until end of cdata + while (text[0] != Ch(']') || text[1] != Ch(']') || text[2] != Ch('>')) + { + if (!text[0]) + RAPIDXML_PARSE_ERROR("unexpected end of data", text); + ++text; + } + text += 3; // Skip ]]> + return 0; // Do not produce CDATA node + } + + // Skip until end of cdata + Ch *value = text; + while (text[0] != Ch(']') || text[1] != Ch(']') || text[2] != Ch('>')) + { + if (!text[0]) + RAPIDXML_PARSE_ERROR("unexpected end of data", text); + ++text; + } + + // Create new cdata node + xml_node *cdata = this->allocate_node(node_cdata); + cdata->value(value, text - value); + + // Place zero terminator after value + if (!(Flags & parse_no_string_terminators)) + *text = Ch('\0'); + + text += 3; // Skip ]]> + return cdata; + } + + // Parse element node + template + xml_node *parse_element(Ch *&text) + { + // Create element node + xml_node *element = this->allocate_node(node_element); + + // Extract element name + Ch *name = text; + skip(text); + if (text == name) + RAPIDXML_PARSE_ERROR("expected element name", text); + element->name(name, text - name); + + // Skip whitespace between element name and attributes or > + skip(text); + + // Parse attributes, if any + parse_node_attributes(text, element); + + // Determine ending type + if (*text == Ch('>')) + { + ++text; + parse_node_contents(text, element); + } + else if (*text == Ch('/')) + { + ++text; + if (*text != Ch('>')) + RAPIDXML_PARSE_ERROR("expected >", text); + ++text; + } + else + RAPIDXML_PARSE_ERROR("expected >", text); + + // Place zero terminator after name + if (!(Flags & parse_no_string_terminators)) + element->name()[element->name_size()] = Ch('\0'); + + // Return parsed element + return element; + } + + // Determine node type, and parse it + template + xml_node *parse_node(Ch *&text) + { + // Parse proper node type + switch (text[0]) + { + + // <... + default: + // Parse and append element node + return parse_element(text); + + // (text); + } + else + { + // Parse PI + return parse_pi(text); + } + + // (text); + } + break; + + // (text); + } + break; + + // (text); + } + + } // switch + + // Attempt to skip other, unrecognized node types starting with ')) + { + if (*text == 0) + RAPIDXML_PARSE_ERROR("unexpected end of data", text); + ++text; + } + ++text; // Skip '>' + return 0; // No node recognized + + } + } + + // Parse contents of the node - children, data etc. + template + void parse_node_contents(Ch *&text, xml_node *node) + { + // For all children and text + while (1) + { + // Skip whitespace between > and node contents + Ch *contents_start = text; // Store start of node contents before whitespace is skipped + skip(text); + Ch next_char = *text; + + // After data nodes, instead of continuing the loop, control jumps here. + // This is because zero termination inside parse_and_append_data() function + // would wreak havoc with the above code. + // Also, skipping whitespace after data nodes is unnecessary. + after_data_node: + + // Determine what comes next: node closing, child node, data node, or 0? + switch (next_char) + { + + // Node closing or child node + case Ch('<'): + if (text[1] == Ch('/')) + { + // Node closing + text += 2; // Skip '(text); + if (!internal::compare(node->name(), node->name_size(), closing_name, text - closing_name, true)) + RAPIDXML_PARSE_ERROR("invalid closing tag name", text); + } + else + { + // No validation, just skip name + skip(text); + } + // Skip remaining whitespace after node name + skip(text); + if (*text != Ch('>')) + RAPIDXML_PARSE_ERROR("expected >", text); + ++text; // Skip '>' + return; // Node closed, finished parsing contents + } + else + { + // Child node + ++text; // Skip '<' + if (xml_node *child = parse_node(text)) + node->append_node(child); + } + break; + + // End of data - error + case Ch('\0'): + RAPIDXML_PARSE_ERROR("unexpected end of data", text); + + // Data node + default: + next_char = parse_and_append_data(node, text, contents_start); + goto after_data_node; // Bypass regular processing after data nodes + + } + } + } + + // Parse XML attributes of the node + template + void parse_node_attributes(Ch *&text, xml_node *node) + { + // For all attributes + while (attribute_name_pred::test(*text)) + { + // Extract attribute name + Ch *name = text; + ++text; // Skip first character of attribute name + skip(text); + if (text == name) + RAPIDXML_PARSE_ERROR("expected attribute name", name); + + // Create new attribute + xml_attribute *attribute = this->allocate_attribute(); + attribute->name(name, text - name); + node->append_attribute(attribute); + + // Skip whitespace after attribute name + skip(text); + + // Skip = + if (*text != Ch('=')) + RAPIDXML_PARSE_ERROR("expected =", text); + ++text; + + // Add terminating zero after name + if (!(Flags & parse_no_string_terminators)) + attribute->name()[attribute->name_size()] = 0; + + // Skip whitespace after = + skip(text); + + // Skip quote and remember if it was ' or " + Ch quote = *text; + if (quote != Ch('\'') && quote != Ch('"')) + RAPIDXML_PARSE_ERROR("expected ' or \"", text); + ++text; + + // Extract attribute value and expand char refs in it + Ch *value = text, *end; + const int AttFlags = Flags & ~parse_normalize_whitespace; // No whitespace normalization in attributes + if (quote == Ch('\'')) + end = skip_and_expand_character_refs, attribute_value_pure_pred, AttFlags>(text); + else + end = skip_and_expand_character_refs, attribute_value_pure_pred, AttFlags>(text); + + // Set attribute value + attribute->value(value, end - value); + + // Make sure that end quote is present + if (*text != quote) + RAPIDXML_PARSE_ERROR("expected ' or \"", text); + ++text; // Skip quote + + // Add terminating zero after value + if (!(Flags & parse_no_string_terminators)) + attribute->value()[attribute->value_size()] = 0; + + // Skip whitespace after attribute value + skip(text); + } + } + + }; + + //! \cond internal + namespace internal + { + + // Whitespace (space \n \r \t) + template + const unsigned char lookup_tables::lookup_whitespace[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, // 0 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 1 + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 2 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 3 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 4 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 5 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 6 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 7 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 8 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // 9 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // A + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // B + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // C + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // D + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // E + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 // F + }; + + // Node name (anything but space \n \r \t / > ? \0) + template + const unsigned char lookup_tables::lookup_node_name[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, // 0 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1 + 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, // 2 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, // 3 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 4 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 5 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 6 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 7 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 8 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 9 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // A + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // B + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // C + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // D + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // E + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 // F + }; + + // Text (i.e. PCDATA) (anything but < \0) + template + const unsigned char lookup_tables::lookup_text[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 0 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 2 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, // 3 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 4 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 5 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 6 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 7 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 8 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 9 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // A + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // B + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // C + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // D + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // E + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 // F + }; + + // Text (i.e. PCDATA) that does not require processing when ws normalization is disabled + // (anything but < \0 &) + template + const unsigned char lookup_tables::lookup_text_pure_no_ws[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 0 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1 + 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 2 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, // 3 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 4 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 5 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 6 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 7 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 8 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 9 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // A + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // B + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // C + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // D + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // E + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 // F + }; + + // Text (i.e. PCDATA) that does not require processing when ws normalizationis is enabled + // (anything but < \0 & space \n \r \t) + template + const unsigned char lookup_tables::lookup_text_pure_with_ws[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, // 0 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1 + 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 2 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, // 3 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 4 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 5 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 6 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 7 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 8 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 9 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // A + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // B + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // C + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // D + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // E + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 // F + }; + + // Attribute name (anything but space \n \r \t / < > = ? ! \0) + template + const unsigned char lookup_tables::lookup_attribute_name[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, // 0 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1 + 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, // 2 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, // 3 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 4 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 5 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 6 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 7 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 8 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 9 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // A + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // B + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // C + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // D + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // E + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 // F + }; + + // Attribute data with single quote (anything but ' \0) + template + const unsigned char lookup_tables::lookup_attribute_data_1[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 0 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1 + 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, // 2 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 3 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 4 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 5 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 6 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 7 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 8 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 9 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // A + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // B + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // C + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // D + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // E + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 // F + }; + + // Attribute data with single quote that does not require processing (anything but ' \0 &) + template + const unsigned char lookup_tables::lookup_attribute_data_1_pure[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 0 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1 + 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, // 2 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 3 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 4 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 5 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 6 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 7 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 8 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 9 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // A + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // B + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // C + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // D + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // E + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 // F + }; + + // Attribute data with double quote (anything but " \0) + template + const unsigned char lookup_tables::lookup_attribute_data_2[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 0 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1 + 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 2 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 3 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 4 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 5 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 6 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 7 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 8 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 9 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // A + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // B + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // C + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // D + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // E + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 // F + }; + + // Attribute data with double quote that does not require processing (anything but " \0 &) + template + const unsigned char lookup_tables::lookup_attribute_data_2_pure[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 0 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 1 + 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 2 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 3 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 4 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 5 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 6 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 7 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 8 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // 9 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // A + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // B + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // C + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // D + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, // E + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 // F + }; + + // Digits (dec and hex, 255 denotes end of numeric character reference) + template + const unsigned char lookup_tables::lookup_digits[256] = + { + // 0 1 2 3 4 5 6 7 8 9 A B C D E F + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // 0 + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // 1 + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // 2 + 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,255,255,255,255,255,255, // 3 + 255, 10, 11, 12, 13, 14, 15,255,255,255,255,255,255,255,255,255, // 4 + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // 5 + 255, 10, 11, 12, 13, 14, 15,255,255,255,255,255,255,255,255,255, // 6 + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // 7 + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // 8 + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // 9 + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // A + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // B + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // C + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // D + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255, // E + 255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255 // F + }; + + // Upper case conversion + template + const unsigned char lookup_tables::lookup_upcase[256] = + { + // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A B C D E F + 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, // 0 + 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, // 1 + 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, // 2 + 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, // 3 + 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, // 4 + 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, // 5 + 96, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, // 6 + 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 123,124,125,126,127, // 7 + 128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143, // 8 + 144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159, // 9 + 160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175, // A + 176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191, // B + 192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207, // C + 208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223, // D + 224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239, // E + 240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255 // F + }; + } + //! \endcond + +} + +// Undefine internal macros +#undef RAPIDXML_PARSE_ERROR + +// On MSVC, restore warnings state +#ifdef _MSC_VER + #pragma warning(pop) +#endif + +#endif diff --git a/lib/include/rapidxml/rapidxml_print.hpp b/lib/include/rapidxml/rapidxml_print.hpp new file mode 100644 index 0000000..4432273 --- /dev/null +++ b/lib/include/rapidxml/rapidxml_print.hpp @@ -0,0 +1,424 @@ +#ifndef RAPIDXML_PRINT_HPP_INCLUDED +#define RAPIDXML_PRINT_HPP_INCLUDED + +// Copyright (C) 2006, 2009 Marcin Kalicinski +// Version 1.13 +// Revision $DateTime: 2009/05/13 01:46:17 $ +//! \file rapidxml_print.hpp This file contains rapidxml printer implementation + +#include "rapidxml.hpp" + +// Only include streams if not disabled +#ifndef RAPIDXML_NO_STREAMS + #include + #include +#endif + +namespace rapidxml +{ + + /////////////////////////////////////////////////////////////////////// + // Printing flags + + const int print_no_indenting = 0x1; //!< Printer flag instructing the printer to suppress indenting of XML. See print() function. + + /////////////////////////////////////////////////////////////////////// + // Internal + + //! \cond internal + namespace internal + { + + /////////////////////////////////////////////////////////////////////////// + // Internal character operations + + // Copy characters from given range to given output iterator + template + inline OutIt copy_chars(const Ch *begin, const Ch *end, OutIt out) + { + while (begin != end) + *out++ = *begin++; + return out; + } + + // Copy characters from given range to given output iterator and expand + // characters into references (< > ' " &) + template + inline OutIt copy_and_expand_chars(const Ch *begin, const Ch *end, Ch noexpand, OutIt out) + { + while (begin != end) + { + if (*begin == noexpand) + { + *out++ = *begin; // No expansion, copy character + } + else + { + switch (*begin) + { + case Ch('<'): + *out++ = Ch('&'); *out++ = Ch('l'); *out++ = Ch('t'); *out++ = Ch(';'); + break; + case Ch('>'): + *out++ = Ch('&'); *out++ = Ch('g'); *out++ = Ch('t'); *out++ = Ch(';'); + break; + case Ch('\''): + *out++ = Ch('&'); *out++ = Ch('a'); *out++ = Ch('p'); *out++ = Ch('o'); *out++ = Ch('s'); *out++ = Ch(';'); + break; + case Ch('"'): + *out++ = Ch('&'); *out++ = Ch('q'); *out++ = Ch('u'); *out++ = Ch('o'); *out++ = Ch('t'); *out++ = Ch(';'); + break; + case Ch('&'): + *out++ = Ch('&'); *out++ = Ch('a'); *out++ = Ch('m'); *out++ = Ch('p'); *out++ = Ch(';'); + break; + default: + *out++ = *begin; // No expansion, copy character + } + } + ++begin; // Step to next character + } + return out; + } + + // Fill given output iterator with repetitions of the same character + template + inline OutIt fill_chars(OutIt out, int n, Ch ch) + { + for (int i = 0; i < n; ++i) + *out++ = ch; + return out; + } + + // Find character + template + inline bool find_char(const Ch *begin, const Ch *end) + { + while (begin != end) + if (*begin++ == ch) + return true; + return false; + } + + /////////////////////////////////////////////////////////////////////////// + // Internal printing operations + + template + inline OutIt print_node(OutIt out, const xml_node *node, int flags, int indent); + + // Print children of the node + template + inline OutIt print_children(OutIt out, const xml_node *node, int flags, int indent) + { + for (xml_node *child = node->first_node(); child; child = child->next_sibling()) + out = print_node(out, child, flags, indent); + return out; + } + + // Print attributes of the node + template + inline OutIt print_attributes(OutIt out, const xml_node *node, int flags) + { + for (xml_attribute *attribute = node->first_attribute(); attribute; attribute = attribute->next_attribute()) + { + if (attribute->name() && attribute->value()) + { + // Print attribute name + *out = Ch(' '), ++out; + out = copy_chars(attribute->name(), attribute->name() + attribute->name_size(), out); + *out = Ch('='), ++out; + // Print attribute value using appropriate quote type + if (find_char(attribute->value(), attribute->value() + attribute->value_size())) + { + *out = Ch('\''), ++out; + out = copy_and_expand_chars(attribute->value(), attribute->value() + attribute->value_size(), Ch('"'), out); + *out = Ch('\''), ++out; + } + else + { + *out = Ch('"'), ++out; + out = copy_and_expand_chars(attribute->value(), attribute->value() + attribute->value_size(), Ch('\''), out); + *out = Ch('"'), ++out; + } + } + } + return out; + } + + // Print data node + template + inline OutIt print_data_node(OutIt out, const xml_node *node, int flags, int indent) + { + assert(node->type() == node_data); + if (!(flags & print_no_indenting)) + out = fill_chars(out, indent, Ch('\t')); + out = copy_and_expand_chars(node->value(), node->value() + node->value_size(), Ch(0), out); + return out; + } + + // Print data node + template + inline OutIt print_cdata_node(OutIt out, const xml_node *node, int flags, int indent) + { + assert(node->type() == node_cdata); + if (!(flags & print_no_indenting)) + out = fill_chars(out, indent, Ch('\t')); + *out = Ch('<'); ++out; + *out = Ch('!'); ++out; + *out = Ch('['); ++out; + *out = Ch('C'); ++out; + *out = Ch('D'); ++out; + *out = Ch('A'); ++out; + *out = Ch('T'); ++out; + *out = Ch('A'); ++out; + *out = Ch('['); ++out; + out = copy_chars(node->value(), node->value() + node->value_size(), out); + *out = Ch(']'); ++out; + *out = Ch(']'); ++out; + *out = Ch('>'); ++out; + return out; + } + + // Print element node + template + inline OutIt print_element_node(OutIt out, const xml_node *node, int flags, int indent) + { + assert(node->type() == node_element); + + // Print element name and attributes, if any + if (!(flags & print_no_indenting)) + out = fill_chars(out, indent, Ch('\t')); + *out = Ch('<'), ++out; + out = copy_chars(node->name(), node->name() + node->name_size(), out); + out = print_attributes(out, node, flags); + + // If node is childless + if (node->value_size() == 0 && !node->first_node()) + { + // Print childless node tag ending + *out = Ch('/'), ++out; + *out = Ch('>'), ++out; + } + else + { + // Print normal node tag ending + *out = Ch('>'), ++out; + + // Test if node contains a single data node only (and no other nodes) + xml_node *child = node->first_node(); + if (!child) + { + // If node has no children, only print its value without indenting + out = copy_and_expand_chars(node->value(), node->value() + node->value_size(), Ch(0), out); + } + else if (child->next_sibling() == 0 && child->type() == node_data) + { + // If node has a sole data child, only print its value without indenting + out = copy_and_expand_chars(child->value(), child->value() + child->value_size(), Ch(0), out); + } + else + { + // Print all children with full indenting + if (!(flags & print_no_indenting)) + *out = Ch('\n'), ++out; + out = print_children(out, node, flags, indent + 1); + if (!(flags & print_no_indenting)) + out = fill_chars(out, indent, Ch('\t')); + } + + // Print node end + *out = Ch('<'), ++out; + *out = Ch('/'), ++out; + out = copy_chars(node->name(), node->name() + node->name_size(), out); + *out = Ch('>'), ++out; + } + return out; + } + + // Print declaration node + template + inline OutIt print_declaration_node(OutIt out, const xml_node *node, int flags, int indent) + { + // Print declaration start + if (!(flags & print_no_indenting)) + out = fill_chars(out, indent, Ch('\t')); + *out = Ch('<'), ++out; + *out = Ch('?'), ++out; + *out = Ch('x'), ++out; + *out = Ch('m'), ++out; + *out = Ch('l'), ++out; + + // Print attributes + out = print_attributes(out, node, flags); + + // Print declaration end + *out = Ch('?'), ++out; + *out = Ch('>'), ++out; + + return out; + } + + // Print comment node + template + inline OutIt print_comment_node(OutIt out, const xml_node *node, int flags, int indent) + { + assert(node->type() == node_comment); + if (!(flags & print_no_indenting)) + out = fill_chars(out, indent, Ch('\t')); + *out = Ch('<'), ++out; + *out = Ch('!'), ++out; + *out = Ch('-'), ++out; + *out = Ch('-'), ++out; + out = copy_chars(node->value(), node->value() + node->value_size(), out); + *out = Ch('-'), ++out; + *out = Ch('-'), ++out; + *out = Ch('>'), ++out; + return out; + } + + // Print doctype node + template + inline OutIt print_doctype_node(OutIt out, const xml_node *node, int flags, int indent) + { + assert(node->type() == node_doctype); + if (!(flags & print_no_indenting)) + out = fill_chars(out, indent, Ch('\t')); + *out = Ch('<'), ++out; + *out = Ch('!'), ++out; + *out = Ch('D'), ++out; + *out = Ch('O'), ++out; + *out = Ch('C'), ++out; + *out = Ch('T'), ++out; + *out = Ch('Y'), ++out; + *out = Ch('P'), ++out; + *out = Ch('E'), ++out; + *out = Ch(' '), ++out; + out = copy_chars(node->value(), node->value() + node->value_size(), out); + *out = Ch('>'), ++out; + return out; + } + + // Print pi node + template + inline OutIt print_pi_node(OutIt out, const xml_node *node, int flags, int indent) + { + assert(node->type() == node_pi); + if (!(flags & print_no_indenting)) + out = fill_chars(out, indent, Ch('\t')); + *out = Ch('<'), ++out; + *out = Ch('?'), ++out; + out = copy_chars(node->name(), node->name() + node->name_size(), out); + *out = Ch(' '), ++out; + out = copy_chars(node->value(), node->value() + node->value_size(), out); + *out = Ch('?'), ++out; + *out = Ch('>'), ++out; + return out; + } + + // Print node + template + inline OutIt print_node(OutIt out, const xml_node *node, int flags, int indent) + { + // Print proper node type + switch (node->type()) + { + + // Document + case node_document: + out = print_children(out, node, flags, indent); + break; + + // Element + case node_element: + out = print_element_node(out, node, flags, indent); + break; + + // Data + case node_data: + out = print_data_node(out, node, flags, indent); + break; + + // CDATA + case node_cdata: + out = print_cdata_node(out, node, flags, indent); + break; + + // Declaration + case node_declaration: + out = print_declaration_node(out, node, flags, indent); + break; + + // Comment + case node_comment: + out = print_comment_node(out, node, flags, indent); + break; + + // Doctype + case node_doctype: + out = print_doctype_node(out, node, flags, indent); + break; + + // Pi + case node_pi: + out = print_pi_node(out, node, flags, indent); + break; + + // Unknown + default: + assert(0); + break; + } + + // If indenting not disabled, add line break after node + if (!(flags & print_no_indenting)) + *out = Ch('\n'), ++out; + + // Return modified iterator + return out; + } + + } + //! \endcond + + /////////////////////////////////////////////////////////////////////////// + // Printing + + //! Prints XML to given output iterator. + //! \param out Output iterator to print to. + //! \param node Node to be printed. Pass xml_document to print entire document. + //! \param flags Flags controlling how XML is printed. + //! \return Output iterator pointing to position immediately after last character of printed text. + template + inline OutIt print(OutIt out, const xml_node &node, int flags = 0) + { + return internal::print_node(out, &node, flags, 0); + } + +#ifndef RAPIDXML_NO_STREAMS + + //! Prints XML to given output stream. + //! \param out Output stream to print to. + //! \param node Node to be printed. Pass xml_document to print entire document. + //! \param flags Flags controlling how XML is printed. + //! \return Output stream. + template + inline std::basic_ostream &print(std::basic_ostream &out, const xml_node &node, int flags = 0) + { + print(std::ostream_iterator(out), node, flags); + return out; + } + + //! Prints formatted XML to given output stream. Uses default printing flags. Use print() function to customize printing process. + //! \param out Output stream to print to. + //! \param node Node to be printed. + //! \return Output stream. + template + inline std::basic_ostream &operator <<(std::basic_ostream &out, const xml_node &node) + { + return print(out, node); + } + +#endif + +} + +#endif diff --git a/lib/include/tnt/jama_cholesky.h b/lib/include/tnt/jama_cholesky.h new file mode 100644 index 0000000..371d2f7 --- /dev/null +++ b/lib/include/tnt/jama_cholesky.h @@ -0,0 +1,258 @@ +#ifndef JAMA_CHOLESKY_H +#define JAMA_CHOLESKY_H + +#include "math.h" + /* needed for sqrt() below. */ + + +namespace JAMA +{ + +using namespace TNT; + +/** +

+ For a symmetric, positive definite matrix A, this function + computes the Cholesky factorization, i.e. it computes a lower + triangular matrix L such that A = L*L'. + If the matrix is not symmetric or positive definite, the function + computes only a partial decomposition. This can be tested with + the is_spd() flag. + +

Typical usage looks like: + 

+	Array2D A(n,n);
+	Array2D L;
+
+	 ... 
+
+	Cholesky chol(A);
+
+	if (chol.is_spd())
+		L = chol.getL();
+		
+  	else
+		cout << "factorization was not complete.\n";
+
+
+ + +

+ (Adapted from JAMA, a Java Matrix Library, developed by jointly + by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). + + */ + +template +class Cholesky +{ + Array2D L_; // lower triangular factor + int isspd; // 1 if matrix to be factored was SPD + +public: + + Cholesky(); + Cholesky(const Array2D &A); + Array2D getL() const; + Array1D solve(const Array1D &B); + Array2D solve(const Array2D &B); + int is_spd() const; + +}; + +template +Cholesky::Cholesky() : L_(0,0), isspd(0) {} + +/** + @return 1, if original matrix to be factored was symmetric + positive-definite (SPD). +*/ +template +int Cholesky::is_spd() const +{ + return isspd; +} + +/** + @return the lower triangular factor, L, such that L*L'=A. +*/ +template +Array2D Cholesky::getL() const +{ + return L_; +} + +/** + Constructs a lower triangular matrix L, such that L*L'= A. + If A is not symmetric positive-definite (SPD), only a + partial factorization is performed. If is_spd() + evalutate true (1) then the factorizaiton was successful. +*/ +template +Cholesky::Cholesky(const Array2D &A) +{ + + + int m = A.dim1(); + int n = A.dim2(); + + isspd = (m == n); + + if (m != n) + { + L_ = Array2D(0,0); + return; + } + + L_ = Array2D(n,n); + + + // Main loop. + for (int j = 0; j < n; j++) + { + Real d(0.0); + for (int k = 0; k < j; k++) + { + Real s(0.0); + for (int i = 0; i < k; i++) + { + s += L_[k][i]*L_[j][i]; + } + L_[j][k] = s = (A[j][k] - s)/L_[k][k]; + d = d + s*s; + isspd = isspd && (A[k][j] == A[j][k]); + } + d = A[j][j] - d; + isspd = isspd && (d > 0.0); + L_[j][j] = sqrt(d > 0.0 ? d : 0.0); + for (int k = j+1; k < n; k++) + { + L_[j][k] = 0.0; + } + } +} + +/** + + Solve a linear system A*x = b, using the previously computed + cholesky factorization of A: L*L'. + + @param B A Matrix with as many rows as A and any number of columns. + @return x so that L*L'*x = b. If b is nonconformat, or if A + was not symmetric posidtive definite, a null (0x0) + array is returned. +*/ +template +Array1D Cholesky::solve(const Array1D &b) +{ + int n = L_.dim1(); + if (b.dim1() != n) + return Array1D(); + + + Array1D x = b.copy(); + + + // Solve L*y = b; + for (int k = 0; k < n; k++) + { + for (int i = 0; i < k; i++) + x[k] -= x[i]*L_[k][i]; + x[k] /= L_[k][k]; + + } + + // Solve L'*X = Y; + for (int k = n-1; k >= 0; k--) + { + for (int i = k+1; i < n; i++) + x[k] -= x[i]*L_[i][k]; + x[k] /= L_[k][k]; + } + + return x; +} + + +/** + + Solve a linear system A*X = B, using the previously computed + cholesky factorization of A: L*L'. + + @param B A Matrix with as many rows as A and any number of columns. + @return X so that L*L'*X = B. If B is nonconformat, or if A + was not symmetric posidtive definite, a null (0x0) + array is returned. +*/ +template +Array2D Cholesky::solve(const Array2D &B) +{ + int n = L_.dim1(); + if (B.dim1() != n) + return Array2D(); + + + Array2D X = B.copy(); + int nx = B.dim2(); + +// Cleve's original code +#if 0 + // Solve L*Y = B; + for (int k = 0; k < n; k++) { + for (int i = k+1; i < n; i++) { + for (int j = 0; j < nx; j++) { + X[i][j] -= X[k][j]*L_[k][i]; + } + } + for (int j = 0; j < nx; j++) { + X[k][j] /= L_[k][k]; + } + } + + // Solve L'*X = Y; + for (int k = n-1; k >= 0; k--) { + for (int j = 0; j < nx; j++) { + X[k][j] /= L_[k][k]; + } + for (int i = 0; i < k; i++) { + for (int j = 0; j < nx; j++) { + X[i][j] -= X[k][j]*L_[k][i]; + } + } + } +#endif + + + // Solve L*y = b; + for (int j=0; j< nx; j++) + { + for (int k = 0; k < n; k++) + { + for (int i = 0; i < k; i++) + X[k][j] -= X[i][j]*L_[k][i]; + X[k][j] /= L_[k][k]; + } + } + + // Solve L'*X = Y; + for (int j=0; j= 0; k--) + { + for (int i = k+1; i < n; i++) + X[k][j] -= X[i][j]*L_[i][k]; + X[k][j] /= L_[k][k]; + } + } + + + + return X; +} + + +} +// namespace JAMA + +#endif +// JAMA_CHOLESKY_H diff --git a/lib/include/tnt/jama_eig.h b/lib/include/tnt/jama_eig.h new file mode 100644 index 0000000..8e3572f --- /dev/null +++ b/lib/include/tnt/jama_eig.h @@ -0,0 +1,1034 @@ +#ifndef JAMA_EIG_H +#define JAMA_EIG_H + + +#include "tnt_array1d.h" +#include "tnt_array2d.h" +#include "tnt_math_utils.h" + +#include +// for min(), max() below + +#include +// for abs() below + +using namespace TNT; +using namespace std; + +// Modification by Willem Jan Palenstijn, 2010-03-11: +// Use std::min() instead of min(), std::max() instead of max() + + +namespace JAMA +{ + +/** + + Computes eigenvalues and eigenvectors of a real (non-complex) + matrix. +

+ If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is + diagonal and the eigenvector matrix V is orthogonal. That is, + the diagonal values of D are the eigenvalues, and + V*V' = I, where I is the identity matrix. The columns of V + represent the eigenvectors in the sense that A*V = V*D. + +

+ If A is not symmetric, then the eigenvalue matrix D is block diagonal + with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, + a + i*b, in 2-by-2 blocks, [a, b; -b, a]. That is, if the complex + eigenvalues look like +

+
+          u + iv     .        .          .      .    .
+            .      u - iv     .          .      .    .
+            .        .      a + ib       .      .    .
+            .        .        .        a - ib   .    .
+            .        .        .          .      x    .
+            .        .        .          .      .    y
+
+ then D looks like +
+
+            u        v        .          .      .    .
+           -v        u        .          .      .    . 
+            .        .        a          b      .    .
+            .        .       -b          a      .    .
+            .        .        .          .      x    .
+            .        .        .          .      .    y
+
+ This keeps V a real matrix in both symmetric and non-symmetric + cases, and A*V = V*D. + + + +

+ The matrix V may be badly + conditioned, or even singular, so the validity of the equation + A = V*D*inverse(V) depends upon the condition number of V. + +

+ (Adapted from JAMA, a Java Matrix Library, developed by jointly + by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). +**/ + +template +class Eigenvalue +{ + + + /** Row and column dimension (square matrix). */ + int n; + + int issymmetric; /* boolean*/ + + /** Arrays for internal storage of eigenvalues. */ + + TNT::Array1D d; /* real part */ + TNT::Array1D e; /* img part */ + + /** Array for internal storage of eigenvectors. */ + TNT::Array2D V; + + /** Array for internal storage of nonsymmetric Hessenberg form. + @serial internal storage of nonsymmetric Hessenberg form. + */ + TNT::Array2D H; + + + /** Working storage for nonsymmetric algorithm. + @serial working storage for nonsymmetric algorithm. + */ + TNT::Array1D ort; + + + // Symmetric Householder reduction to tridiagonal form. + + void tred2() { + + // This is derived from the Algol procedures tred2 by + // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for + // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutine in EISPACK. + + for (int j = 0; j < n; j++) { + d[j] = V[n-1][j]; + } + + // Householder reduction to tridiagonal form. + + for (int i = n-1; i > 0; i--) { + + // Scale to avoid under/overflow. + + Real scale = 0.0; + Real h = 0.0; + for (int k = 0; k < i; k++) { + scale = scale + abs(d[k]); + } + if (scale == 0.0) { + e[i] = d[i-1]; + for (int j = 0; j < i; j++) { + d[j] = V[i-1][j]; + V[i][j] = 0.0; + V[j][i] = 0.0; + } + } else { + + // Generate Householder vector. + + for (int k = 0; k < i; k++) { + d[k] /= scale; + h += d[k] * d[k]; + } + Real f = d[i-1]; + Real g = sqrt(h); + if (f > 0) { + g = -g; + } + e[i] = scale * g; + h = h - f * g; + d[i-1] = f - g; + for (int j = 0; j < i; j++) { + e[j] = 0.0; + } + + // Apply similarity transformation to remaining columns. + + for (int j = 0; j < i; j++) { + f = d[j]; + V[j][i] = f; + g = e[j] + V[j][j] * f; + for (int k = j+1; k <= i-1; k++) { + g += V[k][j] * d[k]; + e[k] += V[k][j] * f; + } + e[j] = g; + } + f = 0.0; + for (int j = 0; j < i; j++) { + e[j] /= h; + f += e[j] * d[j]; + } + Real hh = f / (h + h); + for (int j = 0; j < i; j++) { + e[j] -= hh * d[j]; + } + for (int j = 0; j < i; j++) { + f = d[j]; + g = e[j]; + for (int k = j; k <= i-1; k++) { + V[k][j] -= (f * e[k] + g * d[k]); + } + d[j] = V[i-1][j]; + V[i][j] = 0.0; + } + } + d[i] = h; + } + + // Accumulate transformations. + + for (int i = 0; i < n-1; i++) { + V[n-1][i] = V[i][i]; + V[i][i] = 1.0; + Real h = d[i+1]; + if (h != 0.0) { + for (int k = 0; k <= i; k++) { + d[k] = V[k][i+1] / h; + } + for (int j = 0; j <= i; j++) { + Real g = 0.0; + for (int k = 0; k <= i; k++) { + g += V[k][i+1] * V[k][j]; + } + for (int k = 0; k <= i; k++) { + V[k][j] -= g * d[k]; + } + } + } + for (int k = 0; k <= i; k++) { + V[k][i+1] = 0.0; + } + } + for (int j = 0; j < n; j++) { + d[j] = V[n-1][j]; + V[n-1][j] = 0.0; + } + V[n-1][n-1] = 1.0; + e[0] = 0.0; + } + + // Symmetric tridiagonal QL algorithm. + + void tql2 () { + + // This is derived from the Algol procedures tql2, by + // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for + // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutine in EISPACK. + + for (int i = 1; i < n; i++) { + e[i-1] = e[i]; + } + e[n-1] = 0.0; + + Real f = 0.0; + Real tst1 = 0.0; + Real eps = pow(2.0,-52.0); + for (int l = 0; l < n; l++) { + + // Find small subdiagonal element + + tst1 = std::max(tst1,abs(d[l]) + abs(e[l])); + int m = l; + + // Original while-loop from Java code + while (m < n) { + if (abs(e[m]) <= eps*tst1) { + break; + } + m++; + } + + + // If m == l, d[l] is an eigenvalue, + // otherwise, iterate. + + if (m > l) { + int iter = 0; + do { + iter = iter + 1; // (Could check iteration count here.) + + // Compute implicit shift + + Real g = d[l]; + Real p = (d[l+1] - g) / (2.0 * e[l]); + Real r = hypot(p,1.0); + if (p < 0) { + r = -r; + } + d[l] = e[l] / (p + r); + d[l+1] = e[l] * (p + r); + Real dl1 = d[l+1]; + Real h = g - d[l]; + for (int i = l+2; i < n; i++) { + d[i] -= h; + } + f = f + h; + + // Implicit QL transformation. + + p = d[m]; + Real c = 1.0; + Real c2 = c; + Real c3 = c; + Real el1 = e[l+1]; + Real s = 0.0; + Real s2 = 0.0; + for (int i = m-1; i >= l; i--) { + c3 = c2; + c2 = c; + s2 = s; + g = c * e[i]; + h = c * p; + r = hypot(p,e[i]); + e[i+1] = s * r; + s = e[i] / r; + c = p / r; + p = c * d[i] - s * g; + d[i+1] = h + s * (c * g + s * d[i]); + + // Accumulate transformation. + + for (int k = 0; k < n; k++) { + h = V[k][i+1]; + V[k][i+1] = s * V[k][i] + c * h; + V[k][i] = c * V[k][i] - s * h; + } + } + p = -s * s2 * c3 * el1 * e[l] / dl1; + e[l] = s * p; + d[l] = c * p; + + // Check for convergence. + + } while (abs(e[l]) > eps*tst1); + } + d[l] = d[l] + f; + e[l] = 0.0; + } + + // Sort eigenvalues and corresponding vectors. + + for (int i = 0; i < n-1; i++) { + int k = i; + Real p = d[i]; + for (int j = i+1; j < n; j++) { + if (d[j] < p) { + k = j; + p = d[j]; + } + } + if (k != i) { + d[k] = d[i]; + d[i] = p; + for (int j = 0; j < n; j++) { + p = V[j][i]; + V[j][i] = V[j][k]; + V[j][k] = p; + } + } + } + } + + // Nonsymmetric reduction to Hessenberg form. + + void orthes () { + + // This is derived from the Algol procedures orthes and ortran, + // by Martin and Wilkinson, Handbook for Auto. Comp., + // Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutines in EISPACK. + + int low = 0; + int high = n-1; + + for (int m = low+1; m <= high-1; m++) { + + // Scale column. + + Real scale = 0.0; + for (int i = m; i <= high; i++) { + scale = scale + abs(H[i][m-1]); + } + if (scale != 0.0) { + + // Compute Householder transformation. + + Real h = 0.0; + for (int i = high; i >= m; i--) { + ort[i] = H[i][m-1]/scale; + h += ort[i] * ort[i]; + } + Real g = sqrt(h); + if (ort[m] > 0) { + g = -g; + } + h = h - ort[m] * g; + ort[m] = ort[m] - g; + + // Apply Householder similarity transformation + // H = (I-u*u'/h)*H*(I-u*u')/h) + + for (int j = m; j < n; j++) { + Real f = 0.0; + for (int i = high; i >= m; i--) { + f += ort[i]*H[i][j]; + } + f = f/h; + for (int i = m; i <= high; i++) { + H[i][j] -= f*ort[i]; + } + } + + for (int i = 0; i <= high; i++) { + Real f = 0.0; + for (int j = high; j >= m; j--) { + f += ort[j]*H[i][j]; + } + f = f/h; + for (int j = m; j <= high; j++) { + H[i][j] -= f*ort[j]; + } + } + ort[m] = scale*ort[m]; + H[m][m-1] = scale*g; + } + } + + // Accumulate transformations (Algol's ortran). + + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + V[i][j] = (i == j ? 1.0 : 0.0); + } + } + + for (int m = high-1; m >= low+1; m--) { + if (H[m][m-1] != 0.0) { + for (int i = m+1; i <= high; i++) { + ort[i] = H[i][m-1]; + } + for (int j = m; j <= high; j++) { + Real g = 0.0; + for (int i = m; i <= high; i++) { + g += ort[i] * V[i][j]; + } + // Double division avoids possible underflow + g = (g / ort[m]) / H[m][m-1]; + for (int i = m; i <= high; i++) { + V[i][j] += g * ort[i]; + } + } + } + } + } + + + // Complex scalar division. + + Real cdivr, cdivi; + void cdiv(Real xr, Real xi, Real yr, Real yi) { + Real r,d; + if (abs(yr) > abs(yi)) { + r = yi/yr; + d = yr + r*yi; + cdivr = (xr + r*xi)/d; + cdivi = (xi - r*xr)/d; + } else { + r = yr/yi; + d = yi + r*yr; + cdivr = (r*xr + xi)/d; + cdivi = (r*xi - xr)/d; + } + } + + + // Nonsymmetric reduction from Hessenberg to real Schur form. + + void hqr2 () { + + // This is derived from the Algol procedure hqr2, + // by Martin and Wilkinson, Handbook for Auto. Comp., + // Vol.ii-Linear Algebra, and the corresponding + // Fortran subroutine in EISPACK. + + // Initialize + + int nn = this->n; + int n = nn-1; + int low = 0; + int high = nn-1; + Real eps = pow(2.0,-52.0); + Real exshift = 0.0; + Real p=0,q=0,r=0,s=0,z=0,t,w,x,y; + + // Store roots isolated by balanc and compute matrix norm + + Real norm = 0.0; + for (int i = 0; i < nn; i++) { + if ((i < low) || (i > high)) { + d[i] = H[i][i]; + e[i] = 0.0; + } + for (int j = std::max(i-1,0); j < nn; j++) { + norm = norm + abs(H[i][j]); + } + } + + // Outer loop over eigenvalue index + + int iter = 0; + while (n >= low) { + + // Look for single small sub-diagonal element + + int l = n; + while (l > low) { + s = abs(H[l-1][l-1]) + abs(H[l][l]); + if (s == 0.0) { + s = norm; + } + if (abs(H[l][l-1]) < eps * s) { + break; + } + l--; + } + + // Check for convergence + // One root found + + if (l == n) { + H[n][n] = H[n][n] + exshift; + d[n] = H[n][n]; + e[n] = 0.0; + n--; + iter = 0; + + // Two roots found + + } else if (l == n-1) { + w = H[n][n-1] * H[n-1][n]; + p = (H[n-1][n-1] - H[n][n]) / 2.0; + q = p * p + w; + z = sqrt(abs(q)); + H[n][n] = H[n][n] + exshift; + H[n-1][n-1] = H[n-1][n-1] + exshift; + x = H[n][n]; + + // Real pair + + if (q >= 0) { + if (p >= 0) { + z = p + z; + } else { + z = p - z; + } + d[n-1] = x + z; + d[n] = d[n-1]; + if (z != 0.0) { + d[n] = x - w / z; + } + e[n-1] = 0.0; + e[n] = 0.0; + x = H[n][n-1]; + s = abs(x) + abs(z); + p = x / s; + q = z / s; + r = sqrt(p * p+q * q); + p = p / r; + q = q / r; + + // Row modification + + for (int j = n-1; j < nn; j++) { + z = H[n-1][j]; + H[n-1][j] = q * z + p * H[n][j]; + H[n][j] = q * H[n][j] - p * z; + } + + // Column modification + + for (int i = 0; i <= n; i++) { + z = H[i][n-1]; + H[i][n-1] = q * z + p * H[i][n]; + H[i][n] = q * H[i][n] - p * z; + } + + // Accumulate transformations + + for (int i = low; i <= high; i++) { + z = V[i][n-1]; + V[i][n-1] = q * z + p * V[i][n]; + V[i][n] = q * V[i][n] - p * z; + } + + // Complex pair + + } else { + d[n-1] = x + p; + d[n] = x + p; + e[n-1] = z; + e[n] = -z; + } + n = n - 2; + iter = 0; + + // No convergence yet + + } else { + + // Form shift + + x = H[n][n]; + y = 0.0; + w = 0.0; + if (l < n) { + y = H[n-1][n-1]; + w = H[n][n-1] * H[n-1][n]; + } + + // Wilkinson's original ad hoc shift + + if (iter == 10) { + exshift += x; + for (int i = low; i <= n; i++) { + H[i][i] -= x; + } + s = abs(H[n][n-1]) + abs(H[n-1][n-2]); + x = y = 0.75 * s; + w = -0.4375 * s * s; + } + + // MATLAB's new ad hoc shift + + if (iter == 30) { + s = (y - x) / 2.0; + s = s * s + w; + if (s > 0) { + s = sqrt(s); + if (y < x) { + s = -s; + } + s = x - w / ((y - x) / 2.0 + s); + for (int i = low; i <= n; i++) { + H[i][i] -= s; + } + exshift += s; + x = y = w = 0.964; + } + } + + iter = iter + 1; // (Could check iteration count here.) + + // Look for two consecutive small sub-diagonal elements + + int m = n-2; + while (m >= l) { + z = H[m][m]; + r = x - z; + s = y - z; + p = (r * s - w) / H[m+1][m] + H[m][m+1]; + q = H[m+1][m+1] - z - r - s; + r = H[m+2][m+1]; + s = abs(p) + abs(q) + abs(r); + p = p / s; + q = q / s; + r = r / s; + if (m == l) { + break; + } + if (abs(H[m][m-1]) * (abs(q) + abs(r)) < + eps * (abs(p) * (abs(H[m-1][m-1]) + abs(z) + + abs(H[m+1][m+1])))) { + break; + } + m--; + } + + for (int i = m+2; i <= n; i++) { + H[i][i-2] = 0.0; + if (i > m+2) { + H[i][i-3] = 0.0; + } + } + + // Double QR step involving rows l:n and columns m:n + + for (int k = m; k <= n-1; k++) { + int notlast = (k != n-1); + if (k != m) { + p = H[k][k-1]; + q = H[k+1][k-1]; + r = (notlast ? H[k+2][k-1] : 0.0); + x = abs(p) + abs(q) + abs(r); + if (x != 0.0) { + p = p / x; + q = q / x; + r = r / x; + } + } + if (x == 0.0) { + break; + } + s = sqrt(p * p + q * q + r * r); + if (p < 0) { + s = -s; + } + if (s != 0) { + if (k != m) { + H[k][k-1] = -s * x; + } else if (l != m) { + H[k][k-1] = -H[k][k-1]; + } + p = p + s; + x = p / s; + y = q / s; + z = r / s; + q = q / p; + r = r / p; + + // Row modification + + for (int j = k; j < nn; j++) { + p = H[k][j] + q * H[k+1][j]; + if (notlast) { + p = p + r * H[k+2][j]; + H[k+2][j] = H[k+2][j] - p * z; + } + H[k][j] = H[k][j] - p * x; + H[k+1][j] = H[k+1][j] - p * y; + } + + // Column modification + + for (int i = 0; i <= std::min(n,k+3); i++) { + p = x * H[i][k] + y * H[i][k+1]; + if (notlast) { + p = p + z * H[i][k+2]; + H[i][k+2] = H[i][k+2] - p * r; + } + H[i][k] = H[i][k] - p; + H[i][k+1] = H[i][k+1] - p * q; + } + + // Accumulate transformations + + for (int i = low; i <= high; i++) { + p = x * V[i][k] + y * V[i][k+1]; + if (notlast) { + p = p + z * V[i][k+2]; + V[i][k+2] = V[i][k+2] - p * r; + } + V[i][k] = V[i][k] - p; + V[i][k+1] = V[i][k+1] - p * q; + } + } // (s != 0) + } // k loop + } // check convergence + } // while (n >= low) + + // Backsubstitute to find vectors of upper triangular form + + if (norm == 0.0) { + return; + } + + for (n = nn-1; n >= 0; n--) { + p = d[n]; + q = e[n]; + + // Real vector + + if (q == 0) { + int l = n; + H[n][n] = 1.0; + for (int i = n-1; i >= 0; i--) { + w = H[i][i] - p; + r = 0.0; + for (int j = l; j <= n; j++) { + r = r + H[i][j] * H[j][n]; + } + if (e[i] < 0.0) { + z = w; + s = r; + } else { + l = i; + if (e[i] == 0.0) { + if (w != 0.0) { + H[i][n] = -r / w; + } else { + H[i][n] = -r / (eps * norm); + } + + // Solve real equations + + } else { + x = H[i][i+1]; + y = H[i+1][i]; + q = (d[i] - p) * (d[i] - p) + e[i] * e[i]; + t = (x * s - z * r) / q; + H[i][n] = t; + if (abs(x) > abs(z)) { + H[i+1][n] = (-r - w * t) / x; + } else { + H[i+1][n] = (-s - y * t) / z; + } + } + + // Overflow control + + t = abs(H[i][n]); + if ((eps * t) * t > 1) { + for (int j = i; j <= n; j++) { + H[j][n] = H[j][n] / t; + } + } + } + } + + // Complex vector + + } else if (q < 0) { + int l = n-1; + + // Last vector component imaginary so matrix is triangular + + if (abs(H[n][n-1]) > abs(H[n-1][n])) { + H[n-1][n-1] = q / H[n][n-1]; + H[n-1][n] = -(H[n][n] - p) / H[n][n-1]; + } else { + cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q); + H[n-1][n-1] = cdivr; + H[n-1][n] = cdivi; + } + H[n][n-1] = 0.0; + H[n][n] = 1.0; + for (int i = n-2; i >= 0; i--) { + Real ra,sa,vr,vi; + ra = 0.0; + sa = 0.0; + for (int j = l; j <= n; j++) { + ra = ra + H[i][j] * H[j][n-1]; + sa = sa + H[i][j] * H[j][n]; + } + w = H[i][i] - p; + + if (e[i] < 0.0) { + z = w; + r = ra; + s = sa; + } else { + l = i; + if (e[i] == 0) { + cdiv(-ra,-sa,w,q); + H[i][n-1] = cdivr; + H[i][n] = cdivi; + } else { + + // Solve complex equations + + x = H[i][i+1]; + y = H[i+1][i]; + vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; + vi = (d[i] - p) * 2.0 * q; + if ((vr == 0.0) && (vi == 0.0)) { + vr = eps * norm * (abs(w) + abs(q) + + abs(x) + abs(y) + abs(z)); + } + cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi); + H[i][n-1] = cdivr; + H[i][n] = cdivi; + if (abs(x) > (abs(z) + abs(q))) { + H[i+1][n-1] = (-ra - w * H[i][n-1] + q * H[i][n]) / x; + H[i+1][n] = (-sa - w * H[i][n] - q * H[i][n-1]) / x; + } else { + cdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q); + H[i+1][n-1] = cdivr; + H[i+1][n] = cdivi; + } + } + + // Overflow control + + t = std::max(abs(H[i][n-1]),abs(H[i][n])); + if ((eps * t) * t > 1) { + for (int j = i; j <= n; j++) { + H[j][n-1] = H[j][n-1] / t; + H[j][n] = H[j][n] / t; + } + } + } + } + } + } + + // Vectors of isolated roots + + for (int i = 0; i < nn; i++) { + if (i < low || i > high) { + for (int j = i; j < nn; j++) { + V[i][j] = H[i][j]; + } + } + } + + // Back transformation to get eigenvectors of original matrix + + for (int j = nn-1; j >= low; j--) { + for (int i = low; i <= high; i++) { + z = 0.0; + for (int k = low; k <= std::min(j,high); k++) { + z = z + V[i][k] * H[k][j]; + } + V[i][j] = z; + } + } + } + +public: + + + /** Check for symmetry, then construct the eigenvalue decomposition + @param A Square real (non-complex) matrix + */ + + Eigenvalue(const TNT::Array2D &A) { + n = A.dim2(); + V = Array2D(n,n); + d = Array1D(n); + e = Array1D(n); + + issymmetric = 1; + for (int j = 0; (j < n) && issymmetric; j++) { + for (int i = 0; (i < n) && issymmetric; i++) { + issymmetric = (A[i][j] == A[j][i]); + } + } + + if (issymmetric) { + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + V[i][j] = A[i][j]; + } + } + + // Tridiagonalize. + tred2(); + + // Diagonalize. + tql2(); + + } else { + H = TNT::Array2D(n,n); + ort = TNT::Array1D(n); + + for (int j = 0; j < n; j++) { + for (int i = 0; i < n; i++) { + H[i][j] = A[i][j]; + } + } + + // Reduce to Hessenberg form. + orthes(); + + // Reduce Hessenberg to real Schur form. + hqr2(); + } + } + + + /** Return the eigenvector matrix + @return V + */ + + void getV (TNT::Array2D &V_) { + V_ = V; + return; + } + + /** Return the real parts of the eigenvalues + @return real(diag(D)) + */ + + void getRealEigenvalues (TNT::Array1D &d_) { + d_ = d; + return ; + } + + /** Return the imaginary parts of the eigenvalues + in parameter e_. + + @pararm e_: new matrix with imaginary parts of the eigenvalues. + */ + void getImagEigenvalues (TNT::Array1D &e_) { + e_ = e; + return; + } + + +/** + Computes the block diagonal eigenvalue matrix. + If the original matrix A is not symmetric, then the eigenvalue + matrix D is block diagonal with the real eigenvalues in 1-by-1 + blocks and any complex eigenvalues, + a + i*b, in 2-by-2 blocks, [a, b; -b, a]. That is, if the complex + eigenvalues look like +

+
+          u + iv     .        .          .      .    .
+            .      u - iv     .          .      .    .
+            .        .      a + ib       .      .    .
+            .        .        .        a - ib   .    .
+            .        .        .          .      x    .
+            .        .        .          .      .    y
+
+ then D looks like +
+
+            u        v        .          .      .    .
+           -v        u        .          .      .    . 
+            .        .        a          b      .    .
+            .        .       -b          a      .    .
+            .        .        .          .      x    .
+            .        .        .          .      .    y
+
+ This keeps V a real matrix in both symmetric and non-symmetric + cases, and A*V = V*D. + + @param D: upon return, the matrix is filled with the block diagonal + eigenvalue matrix. + +*/ + void getD (TNT::Array2D &D) { + D = Array2D(n,n); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + D[i][j] = 0.0; + } + D[i][i] = d[i]; + if (e[i] > 0) { + D[i][i+1] = e[i]; + } else if (e[i] < 0) { + D[i][i-1] = e[i]; + } + } + } +}; + +} //namespace JAMA + + +#endif +// JAMA_EIG_H diff --git a/lib/include/tnt/jama_lu.h b/lib/include/tnt/jama_lu.h new file mode 100644 index 0000000..e95b433 --- /dev/null +++ b/lib/include/tnt/jama_lu.h @@ -0,0 +1,323 @@ +#ifndef JAMA_LU_H +#define JAMA_LU_H + +#include "tnt.h" +#include +//for min(), max() below + +using namespace TNT; +using namespace std; + + +// Modification by Willem Jan Palenstijn, 2010-03-11: +// Use std::min() instead of min() + +namespace JAMA +{ + + /** LU Decomposition. +

+ For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n + unit lower triangular matrix L, an n-by-n upper triangular matrix U, + and a permutation vector piv of length m so that A(piv,:) = L*U. + If m < n, then L is m-by-m and U is m-by-n. +

+ The LU decompostion with pivoting always exists, even if the matrix is + singular, so the constructor will never fail. The primary use of the + LU decomposition is in the solution of square systems of simultaneous + linear equations. This will fail if isNonsingular() returns false. + */ +template +class LU +{ + + + + /* Array for internal storage of decomposition. */ + Array2D LU_; + int m, n, pivsign; + Array1D piv; + + + Array2D permute_copy(const Array2D &A, + const Array1D &piv, int j0, int j1) + { + int piv_length = piv.dim(); + + Array2D X(piv_length, j1-j0+1); + + + for (int i = 0; i < piv_length; i++) + for (int j = j0; j <= j1; j++) + X[i][j-j0] = A[piv[i]][j]; + + return X; + } + + Array1D permute_copy(const Array1D &A, + const Array1D &piv) + { + int piv_length = piv.dim(); + if (piv_length != A.dim()) + return Array1D(); + + Array1D x(piv_length); + + + for (int i = 0; i < piv_length; i++) + x[i] = A[piv[i]]; + + return x; + } + + + public : + + /** LU Decomposition + @param A Rectangular matrix + @return LU Decomposition object to access L, U and piv. + */ + + LU (const Array2D &A) : LU_(A.copy()), m(A.dim1()), n(A.dim2()), + piv(A.dim1()) + + { + + // Use a "left-looking", dot-product, Crout/Doolittle algorithm. + + + for (int i = 0; i < m; i++) { + piv[i] = i; + } + pivsign = 1; + Real *LUrowi = 0;; + Array1D LUcolj(m); + + // Outer loop. + + for (int j = 0; j < n; j++) { + + // Make a copy of the j-th column to localize references. + + for (int i = 0; i < m; i++) { + LUcolj[i] = LU_[i][j]; + } + + // Apply previous transformations. + + for (int i = 0; i < m; i++) { + LUrowi = LU_[i]; + + // Most of the time is spent in the following dot product. + + int kmax = std::min(i,j); + double s = 0.0; + for (int k = 0; k < kmax; k++) { + s += LUrowi[k]*LUcolj[k]; + } + + LUrowi[j] = LUcolj[i] -= s; + } + + // Find pivot and exchange if necessary. + + int p = j; + for (int i = j+1; i < m; i++) { + if (abs(LUcolj[i]) > abs(LUcolj[p])) { + p = i; + } + } + if (p != j) { + int k=0; + for (k = 0; k < n; k++) { + double t = LU_[p][k]; + LU_[p][k] = LU_[j][k]; + LU_[j][k] = t; + } + k = piv[p]; + piv[p] = piv[j]; + piv[j] = k; + pivsign = -pivsign; + } + + // Compute multipliers. + + if ((j < m) && (LU_[j][j] != 0.0)) { + for (int i = j+1; i < m; i++) { + LU_[i][j] /= LU_[j][j]; + } + } + } + } + + + /** Is the matrix nonsingular? + @return 1 (true) if upper triangular factor U (and hence A) + is nonsingular, 0 otherwise. + */ + + int isNonsingular () { + for (int j = 0; j < n; j++) { + if (LU_[j][j] == 0) + return 0; + } + return 1; + } + + /** Return lower triangular factor + @return L + */ + + Array2D getL () { + Array2D L_(m,n); + for (int i = 0; i < m; i++) { + for (int j = 0; j < n; j++) { + if (i > j) { + L_[i][j] = LU_[i][j]; + } else if (i == j) { + L_[i][j] = 1.0; + } else { + L_[i][j] = 0.0; + } + } + } + return L_; + } + + /** Return upper triangular factor + @return U portion of LU factorization. + */ + + Array2D getU () { + Array2D U_(n,n); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + if (i <= j) { + U_[i][j] = LU_[i][j]; + } else { + U_[i][j] = 0.0; + } + } + } + return U_; + } + + /** Return pivot permutation vector + @return piv + */ + + Array1D getPivot () { + return piv; + } + + + /** Compute determinant using LU factors. + @return determinant of A, or 0 if A is not square. + */ + + Real det () { + if (m != n) { + return Real(0); + } + Real d = Real(pivsign); + for (int j = 0; j < n; j++) { + d *= LU_[j][j]; + } + return d; + } + + /** Solve A*X = B + @param B A Matrix with as many rows as A and any number of columns. + @return X so that L*U*X = B(piv,:), if B is nonconformant, returns + 0x0 (null) array. + */ + + Array2D solve (const Array2D &B) + { + + /* Dimensions: A is mxn, X is nxk, B is mxk */ + + if (B.dim1() != m) { + return Array2D(0,0); + } + if (!isNonsingular()) { + return Array2D(0,0); + } + + // Copy right hand side with pivoting + int nx = B.dim2(); + + + Array2D X = permute_copy(B, piv, 0, nx-1); + + // Solve L*Y = B(piv,:) + for (int k = 0; k < n; k++) { + for (int i = k+1; i < n; i++) { + for (int j = 0; j < nx; j++) { + X[i][j] -= X[k][j]*LU_[i][k]; + } + } + } + // Solve U*X = Y; + for (int k = n-1; k >= 0; k--) { + for (int j = 0; j < nx; j++) { + X[k][j] /= LU_[k][k]; + } + for (int i = 0; i < k; i++) { + for (int j = 0; j < nx; j++) { + X[i][j] -= X[k][j]*LU_[i][k]; + } + } + } + return X; + } + + + /** Solve A*x = b, where x and b are vectors of length equal + to the number of rows in A. + + @param b a vector (Array1D> of length equal to the first dimension + of A. + @return x a vector (Array1D> so that L*U*x = b(piv), if B is nonconformant, + returns 0x0 (null) array. + */ + + Array1D solve (const Array1D &b) + { + + /* Dimensions: A is mxn, X is nxk, B is mxk */ + + if (b.dim1() != m) { + return Array1D(); + } + if (!isNonsingular()) { + return Array1D(); + } + + + Array1D x = permute_copy(b, piv); + + // Solve L*Y = B(piv) + for (int k = 0; k < n; k++) { + for (int i = k+1; i < n; i++) { + x[i] -= x[k]*LU_[i][k]; + } + } + + // Solve U*X = Y; + for (int k = n-1; k >= 0; k--) { + x[k] /= LU_[k][k]; + for (int i = 0; i < k; i++) + x[i] -= x[k]*LU_[i][k]; + } + + + return x; + } + +}; /* class LU */ + +} /* namespace JAMA */ + +#endif +/* JAMA_LU_H */ diff --git a/lib/include/tnt/jama_qr.h b/lib/include/tnt/jama_qr.h new file mode 100644 index 0000000..98d37f5 --- /dev/null +++ b/lib/include/tnt/jama_qr.h @@ -0,0 +1,326 @@ +#ifndef JAMA_QR_H +#define JAMA_QR_H + +#include "tnt_array1d.h" +#include "tnt_array2d.h" +#include "tnt_math_utils.h" + +namespace JAMA +{ + +/** +

+ Classical QR Decompisition: + for an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n + orthogonal matrix Q and an n-by-n upper triangular matrix R so that + A = Q*R. +

+ The QR decompostion always exists, even if the matrix does not have + full rank, so the constructor will never fail. The primary use of the + QR decomposition is in the least squares solution of nonsquare systems + of simultaneous linear equations. This will fail if isFullRank() + returns 0 (false). + +

+ The Q and R factors can be retrived via the getQ() and getR() + methods. Furthermore, a solve() method is provided to find the + least squares solution of Ax=b using the QR factors. + +

+ (Adapted from JAMA, a Java Matrix Library, developed by jointly + by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). +*/ + +template +class QR { + + + /** Array for internal storage of decomposition. + @serial internal array storage. + */ + + TNT::Array2D QR_; + + /** Row and column dimensions. + @serial column dimension. + @serial row dimension. + */ + int m, n; + + /** Array for internal storage of diagonal of R. + @serial diagonal of R. + */ + TNT::Array1D Rdiag; + + +public: + +/** + Create a QR factorization object for A. + + @param A rectangular (m>=n) matrix. +*/ + QR(const TNT::Array2D &A) /* constructor */ + { + QR_ = A.copy(); + m = A.dim1(); + n = A.dim2(); + Rdiag = TNT::Array1D(n); + int i=0, j=0, k=0; + + // Main loop. + for (k = 0; k < n; k++) { + // Compute 2-norm of k-th column without under/overflow. + Real nrm = 0; + for (i = k; i < m; i++) { + nrm = TNT::hypot(nrm,QR_[i][k]); + } + + if (nrm != 0.0) { + // Form k-th Householder vector. + if (QR_[k][k] < 0) { + nrm = -nrm; + } + for (i = k; i < m; i++) { + QR_[i][k] /= nrm; + } + QR_[k][k] += 1.0; + + // Apply transformation to remaining columns. + for (j = k+1; j < n; j++) { + Real s = 0.0; + for (i = k; i < m; i++) { + s += QR_[i][k]*QR_[i][j]; + } + s = -s/QR_[k][k]; + for (i = k; i < m; i++) { + QR_[i][j] += s*QR_[i][k]; + } + } + } + Rdiag[k] = -nrm; + } + } + + +/** + Flag to denote the matrix is of full rank. + + @return 1 if matrix is full rank, 0 otherwise. +*/ + int isFullRank() const + { + for (int j = 0; j < n; j++) + { + if (Rdiag[j] == 0) + return 0; + } + return 1; + } + + + + + /** + + Retreive the Householder vectors from QR factorization + @returns lower trapezoidal matrix whose columns define the reflections + */ + + TNT::Array2D getHouseholder (void) const + { + TNT::Array2D H(m,n); + + /* note: H is completely filled in by algorithm, so + initializaiton of H is not necessary. + */ + for (int i = 0; i < m; i++) + { + for (int j = 0; j < n; j++) + { + if (i >= j) { + H[i][j] = QR_[i][j]; + } else { + H[i][j] = 0.0; + } + } + } + return H; + } + + + + /** Return the upper triangular factor, R, of the QR factorization + @return R + */ + + TNT::Array2D getR() const + { + TNT::Array2D R(n,n); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + if (i < j) { + R[i][j] = QR_[i][j]; + } else if (i == j) { + R[i][j] = Rdiag[i]; + } else { + R[i][j] = 0.0; + } + } + } + return R; + } + + + + + + /** + Generate and return the (economy-sized) orthogonal factor + @param Q the (ecnomy-sized) orthogonal factor (Q*R=A). + */ + + TNT::Array2D getQ() const + { + int i=0, j=0, k=0; + + TNT::Array2D Q(m,n); + for (k = n-1; k >= 0; k--) { + for (i = 0; i < m; i++) { + Q[i][k] = 0.0; + } + Q[k][k] = 1.0; + for (j = k; j < n; j++) { + if (QR_[k][k] != 0) { + Real s = 0.0; + for (i = k; i < m; i++) { + s += QR_[i][k]*Q[i][j]; + } + s = -s/QR_[k][k]; + for (i = k; i < m; i++) { + Q[i][j] += s*QR_[i][k]; + } + } + } + } + return Q; + } + + + /** Least squares solution of A*x = b + @param B m-length array (vector). + @return x n-length array (vector) that minimizes the two norm of Q*R*X-B. + If B is non-conformant, or if QR.isFullRank() is false, + the routine returns a null (0-length) vector. + */ + + TNT::Array1D solve(const TNT::Array1D &b) const + { + if (b.dim1() != m) /* arrays must be conformant */ + return TNT::Array1D(); + + if ( !isFullRank() ) /* matrix is rank deficient */ + { + return TNT::Array1D(); + } + + TNT::Array1D x = b.copy(); + + // Compute Y = transpose(Q)*b + for (int k = 0; k < n; k++) + { + Real s = 0.0; + for (int i = k; i < m; i++) + { + s += QR_[i][k]*x[i]; + } + s = -s/QR_[k][k]; + for (int i = k; i < m; i++) + { + x[i] += s*QR_[i][k]; + } + } + // Solve R*X = Y; + for (int k = n-1; k >= 0; k--) + { + x[k] /= Rdiag[k]; + for (int i = 0; i < k; i++) { + x[i] -= x[k]*QR_[i][k]; + } + } + + + /* return n x nx portion of X */ + TNT::Array1D x_(n); + for (int i=0; i solve(const TNT::Array2D &B) const + { + if (B.dim1() != m) /* arrays must be conformant */ + return TNT::Array2D(0,0); + + if ( !isFullRank() ) /* matrix is rank deficient */ + { + return TNT::Array2D(0,0); + } + + int nx = B.dim2(); + TNT::Array2D X = B.copy(); + int i=0, j=0, k=0; + + // Compute Y = transpose(Q)*B + for (k = 0; k < n; k++) { + for (j = 0; j < nx; j++) { + Real s = 0.0; + for (i = k; i < m; i++) { + s += QR_[i][k]*X[i][j]; + } + s = -s/QR_[k][k]; + for (i = k; i < m; i++) { + X[i][j] += s*QR_[i][k]; + } + } + } + // Solve R*X = Y; + for (k = n-1; k >= 0; k--) { + for (j = 0; j < nx; j++) { + X[k][j] /= Rdiag[k]; + } + for (i = 0; i < k; i++) { + for (j = 0; j < nx; j++) { + X[i][j] -= X[k][j]*QR_[i][k]; + } + } + } + + + /* return n x nx portion of X */ + TNT::Array2D X_(n,nx); + for (i=0; i +// for min(), max() below +#include +// for abs() below + +using namespace TNT; +using namespace std; + + +// Modification by Willem Jan Palenstijn, 2010-03-11: +// Use std::min() instead of min(), std::max() instead of max() + + +namespace JAMA +{ + /** Singular Value Decomposition. +

+ For an m-by-n matrix A with m >= n, the singular value decomposition is + an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and + an n-by-n orthogonal matrix V so that A = U*S*V'. +

+ The singular values, sigma[k] = S[k][k], are ordered so that + sigma[0] >= sigma[1] >= ... >= sigma[n-1]. +

+ The singular value decompostion always exists, so the constructor will + never fail. The matrix condition number and the effective numerical + rank can be computed from this decomposition. + +

+ (Adapted from JAMA, a Java Matrix Library, developed by jointly + by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama). + */ +template +class SVD +{ + + + Array2D U, V; + Array1D s; + int m, n; + + public: + + + SVD (const Array2D &Arg) { + + + m = Arg.dim1(); + n = Arg.dim2(); + int nu = std::min(m,n); + s = Array1D(std::min(m+1,n)); + U = Array2D(m, nu, Real(0)); + V = Array2D(n,n); + Array1D e(n); + Array1D work(m); + Array2D A(Arg.copy()); + int wantu = 1; /* boolean */ + int wantv = 1; /* boolean */ + int i=0, j=0, k=0; + + // Reduce A to bidiagonal form, storing the diagonal elements + // in s and the super-diagonal elements in e. + + int nct = std::min(m-1,n); + int nrt = std::max(0,std::min(n-2,m)); + for (k = 0; k < std::max(nct,nrt); k++) { + if (k < nct) { + + // Compute the transformation for the k-th column and + // place the k-th diagonal in s[k]. + // Compute 2-norm of k-th column without under/overflow. + s[k] = 0; + for (i = k; i < m; i++) { + s[k] = hypot(s[k],A[i][k]); + } + if (s[k] != 0.0) { + if (A[k][k] < 0.0) { + s[k] = -s[k]; + } + for (i = k; i < m; i++) { + A[i][k] /= s[k]; + } + A[k][k] += 1.0; + } + s[k] = -s[k]; + } + for (j = k+1; j < n; j++) { + if ((k < nct) && (s[k] != 0.0)) { + + // Apply the transformation. + + Real t(0.0); + for (i = k; i < m; i++) { + t += A[i][k]*A[i][j]; + } + t = -t/A[k][k]; + for (i = k; i < m; i++) { + A[i][j] += t*A[i][k]; + } + } + + // Place the k-th row of A into e for the + // subsequent calculation of the row transformation. + + e[j] = A[k][j]; + } + if (wantu & (k < nct)) { + + // Place the transformation in U for subsequent back + // multiplication. + + for (i = k; i < m; i++) { + U[i][k] = A[i][k]; + } + } + if (k < nrt) { + + // Compute the k-th row transformation and place the + // k-th super-diagonal in e[k]. + // Compute 2-norm without under/overflow. + e[k] = 0; + for (i = k+1; i < n; i++) { + e[k] = hypot(e[k],e[i]); + } + if (e[k] != 0.0) { + if (e[k+1] < 0.0) { + e[k] = -e[k]; + } + for (i = k+1; i < n; i++) { + e[i] /= e[k]; + } + e[k+1] += 1.0; + } + e[k] = -e[k]; + if ((k+1 < m) & (e[k] != 0.0)) { + + // Apply the transformation. + + for (i = k+1; i < m; i++) { + work[i] = 0.0; + } + for (j = k+1; j < n; j++) { + for (i = k+1; i < m; i++) { + work[i] += e[j]*A[i][j]; + } + } + for (j = k+1; j < n; j++) { + Real t(-e[j]/e[k+1]); + for (i = k+1; i < m; i++) { + A[i][j] += t*work[i]; + } + } + } + if (wantv) { + + // Place the transformation in V for subsequent + // back multiplication. + + for (i = k+1; i < n; i++) { + V[i][k] = e[i]; + } + } + } + } + + // Set up the final bidiagonal matrix or order p. + + int p = std::min(n,m+1); + if (nct < n) { + s[nct] = A[nct][nct]; + } + if (m < p) { + s[p-1] = 0.0; + } + if (nrt+1 < p) { + e[nrt] = A[nrt][p-1]; + } + e[p-1] = 0.0; + + // If required, generate U. + + if (wantu) { + for (j = nct; j < nu; j++) { + for (i = 0; i < m; i++) { + U[i][j] = 0.0; + } + U[j][j] = 1.0; + } + for (k = nct-1; k >= 0; k--) { + if (s[k] != 0.0) { + for (j = k+1; j < nu; j++) { + Real t(0.0); + for (i = k; i < m; i++) { + t += U[i][k]*U[i][j]; + } + t = -t/U[k][k]; + for (i = k; i < m; i++) { + U[i][j] += t*U[i][k]; + } + } + for (i = k; i < m; i++ ) { + U[i][k] = -U[i][k]; + } + U[k][k] = 1.0 + U[k][k]; + for (i = 0; i < k-1; i++) { + U[i][k] = 0.0; + } + } else { + for (i = 0; i < m; i++) { + U[i][k] = 0.0; + } + U[k][k] = 1.0; + } + } + } + + // If required, generate V. + + if (wantv) { + for (k = n-1; k >= 0; k--) { + if ((k < nrt) & (e[k] != 0.0)) { + for (j = k+1; j < nu; j++) { + Real t(0.0); + for (i = k+1; i < n; i++) { + t += V[i][k]*V[i][j]; + } + t = -t/V[k+1][k]; + for (i = k+1; i < n; i++) { + V[i][j] += t*V[i][k]; + } + } + } + for (i = 0; i < n; i++) { + V[i][k] = 0.0; + } + V[k][k] = 1.0; + } + } + + // Main iteration loop for the singular values. + + int pp = p-1; + int iter = 0; + Real eps(pow(2.0,-52.0)); + while (p > 0) { + int k=0; + int kase=0; + + // Here is where a test for too many iterations would go. + + // This section of the program inspects for + // negligible elements in the s and e arrays. On + // completion the variables kase and k are set as follows. + + // kase = 1 if s(p) and e[k-1] are negligible and k

= -1; k--) { + if (k == -1) { + break; + } + if (abs(e[k]) <= eps*(abs(s[k]) + abs(s[k+1]))) { + e[k] = 0.0; + break; + } + } + if (k == p-2) { + kase = 4; + } else { + int ks; + for (ks = p-1; ks >= k; ks--) { + if (ks == k) { + break; + } + Real t( (ks != p ? abs(e[ks]) : 0.) + + (ks != k+1 ? abs(e[ks-1]) : 0.)); + if (abs(s[ks]) <= eps*t) { + s[ks] = 0.0; + break; + } + } + if (ks == k) { + kase = 3; + } else if (ks == p-1) { + kase = 1; + } else { + kase = 2; + k = ks; + } + } + k++; + + // Perform the task indicated by kase. + + switch (kase) { + + // Deflate negligible s(p). + + case 1: { + Real f(e[p-2]); + e[p-2] = 0.0; + for (j = p-2; j >= k; j--) { + Real t( hypot(s[j],f)); + Real cs(s[j]/t); + Real sn(f/t); + s[j] = t; + if (j != k) { + f = -sn*e[j-1]; + e[j-1] = cs*e[j-1]; + } + if (wantv) { + for (i = 0; i < n; i++) { + t = cs*V[i][j] + sn*V[i][p-1]; + V[i][p-1] = -sn*V[i][j] + cs*V[i][p-1]; + V[i][j] = t; + } + } + } + } + break; + + // Split at negligible s(k). + + case 2: { + Real f(e[k-1]); + e[k-1] = 0.0; + for (j = k; j < p; j++) { + Real t(hypot(s[j],f)); + Real cs( s[j]/t); + Real sn(f/t); + s[j] = t; + f = -sn*e[j]; + e[j] = cs*e[j]; + if (wantu) { + for (i = 0; i < m; i++) { + t = cs*U[i][j] + sn*U[i][k-1]; + U[i][k-1] = -sn*U[i][j] + cs*U[i][k-1]; + U[i][j] = t; + } + } + } + } + break; + + // Perform one qr step. + + case 3: { + + // Calculate the shift. + + Real scale = std::max(std::max(std::max(std::max( + abs(s[p-1]),abs(s[p-2])),abs(e[p-2])), + abs(s[k])),abs(e[k])); + Real sp = s[p-1]/scale; + Real spm1 = s[p-2]/scale; + Real epm1 = e[p-2]/scale; + Real sk = s[k]/scale; + Real ek = e[k]/scale; + Real b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/2.0; + Real c = (sp*epm1)*(sp*epm1); + Real shift = 0.0; + if ((b != 0.0) || (c != 0.0)) { + shift = sqrt(b*b + c); + if (b < 0.0) { + shift = -shift; + } + shift = c/(b + shift); + } + Real f = (sk + sp)*(sk - sp) + shift; + Real g = sk*ek; + + // Chase zeros. + + for (j = k; j < p-1; j++) { + Real t = hypot(f,g); + Real cs = f/t; + Real sn = g/t; + if (j != k) { + e[j-1] = t; + } + f = cs*s[j] + sn*e[j]; + e[j] = cs*e[j] - sn*s[j]; + g = sn*s[j+1]; + s[j+1] = cs*s[j+1]; + if (wantv) { + for (i = 0; i < n; i++) { + t = cs*V[i][j] + sn*V[i][j+1]; + V[i][j+1] = -sn*V[i][j] + cs*V[i][j+1]; + V[i][j] = t; + } + } + t = hypot(f,g); + cs = f/t; + sn = g/t; + s[j] = t; + f = cs*e[j] + sn*s[j+1]; + s[j+1] = -sn*e[j] + cs*s[j+1]; + g = sn*e[j+1]; + e[j+1] = cs*e[j+1]; + if (wantu && (j < m-1)) { + for (i = 0; i < m; i++) { + t = cs*U[i][j] + sn*U[i][j+1]; + U[i][j+1] = -sn*U[i][j] + cs*U[i][j+1]; + U[i][j] = t; + } + } + } + e[p-2] = f; + iter = iter + 1; + } + break; + + // Convergence. + + case 4: { + + // Make the singular values positive. + + if (s[k] <= 0.0) { + s[k] = (s[k] < 0.0 ? -s[k] : 0.0); + if (wantv) { + for (i = 0; i <= pp; i++) { + V[i][k] = -V[i][k]; + } + } + } + + // Order the singular values. + + while (k < pp) { + if (s[k] >= s[k+1]) { + break; + } + Real t = s[k]; + s[k] = s[k+1]; + s[k+1] = t; + if (wantv && (k < n-1)) { + for (i = 0; i < n; i++) { + t = V[i][k+1]; V[i][k+1] = V[i][k]; V[i][k] = t; + } + } + if (wantu && (k < m-1)) { + for (i = 0; i < m; i++) { + t = U[i][k+1]; U[i][k+1] = U[i][k]; U[i][k] = t; + } + } + k++; + } + iter = 0; + p--; + } + break; + } + } + } + + + void getU (Array2D &A) + { + int minm = std::min(m+1,n); + + A = Array2D(m, minm); + + for (int i=0; i &A) + { + A = V; + } + + /** Return the one-dimensional array of singular values */ + + void getSingularValues (Array1D &x) + { + x = s; + } + + /** Return the diagonal matrix of singular values + @return S + */ + + void getS (Array2D &A) { + A = Array2D(n,n); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + A[i][j] = 0.0; + } + A[i][i] = s[i]; + } + } + + /** Two norm (max(S)) */ + + Real norm2 () { + return s[0]; + } + + /** Two norm of condition number (max(S)/min(S)) */ + + Real cond () { + return s[0]/s[std::min(m,n)-1]; + } + + /** Effective numerical matrix rank + @return Number of nonnegligible singular values. + */ + + int rank () + { + Real eps = pow(2.0,-52.0); + Real tol = std::max(m,n)*s[0]*eps; + int r = 0; + for (int i = 0; i < s.dim(); i++) { + if (s[i] > tol) { + r++; + } + } + return r; + } +}; + +} +#endif +// JAMA_SVD_H diff --git a/lib/include/tnt/tnt.h b/lib/include/tnt/tnt.h new file mode 100644 index 0000000..92463e0 --- /dev/null +++ b/lib/include/tnt/tnt.h @@ -0,0 +1,64 @@ +/* +* +* Template Numerical Toolkit (TNT): Linear Algebra Module +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + +#ifndef TNT_H +#define TNT_H + + + +//--------------------------------------------------------------------- +// Define this macro if you want TNT to track some of the out-of-bounds +// indexing. This can encur a small run-time overhead, but is recommended +// while developing code. It can be turned off for production runs. +// +// #define TNT_BOUNDS_CHECK +//--------------------------------------------------------------------- +// + +//#define TNT_BOUNDS_CHECK + + + +#include "tnt_version.h" +#include "tnt_math_utils.h" +#include "tnt_array1d.h" +#include "tnt_array2d.h" +#include "tnt_array3d.h" +#include "tnt_array1d_utils.h" +#include "tnt_array2d_utils.h" +#include "tnt_array3d_utils.h" + +#include "tnt_fortran_array1d.h" +#include "tnt_fortran_array2d.h" +#include "tnt_fortran_array3d.h" +#include "tnt_fortran_array1d_utils.h" +#include "tnt_fortran_array2d_utils.h" +#include "tnt_fortran_array3d_utils.h" + +#include "tnt_sparse_matrix_csr.h" + +#include "tnt_stopwatch.h" +#include "tnt_subscript.h" +#include "tnt_vec.h" +#include "tnt_cmat.h" + + +#endif +// TNT_H diff --git a/lib/include/tnt/tnt_array1d.h b/lib/include/tnt/tnt_array1d.h new file mode 100644 index 0000000..858df57 --- /dev/null +++ b/lib/include/tnt/tnt_array1d.h @@ -0,0 +1,278 @@ +/* +* +* Template Numerical Toolkit (TNT) +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + + +#ifndef TNT_ARRAY1D_H +#define TNT_ARRAY1D_H + +//#include +#include + +#ifdef TNT_BOUNDS_CHECK +#include +#endif + + +#include "tnt_i_refvec.h" + +namespace TNT +{ + +template +class Array1D +{ + + private: + + /* ... */ + i_refvec v_; + int n_; + T* data_; /* this normally points to v_.begin(), but + * could also point to a portion (subvector) + * of v_. + */ + + void copy_(T* p, const T* q, int len) const; + void set_(T* begin, T* end, const T& val); + + + public: + + typedef T value_type; + + + Array1D(); + explicit Array1D(int n); + Array1D(int n, const T &a); + Array1D(int n, T *a); + inline Array1D(const Array1D &A); + inline operator T*(); + inline operator const T*(); + inline Array1D & operator=(const T &a); + inline Array1D & operator=(const Array1D &A); + inline Array1D & ref(const Array1D &A); + Array1D copy() const; + Array1D & inject(const Array1D & A); + inline T& operator[](int i); + inline const T& operator[](int i) const; + inline int dim1() const; + inline int dim() const; + ~Array1D(); + + + /* ... extended interface ... */ + + inline int ref_count() const; + inline Array1D subarray(int i0, int i1); + +}; + + + + +template +Array1D::Array1D() : v_(), n_(0), data_(0) {} + +template +Array1D::Array1D(const Array1D &A) : v_(A.v_), n_(A.n_), + data_(A.data_) +{ +#ifdef TNT_DEBUG + std::cout << "Created Array1D(const Array1D &A) \n"; +#endif + +} + + +template +Array1D::Array1D(int n) : v_(n), n_(n), data_(v_.begin()) +{ +#ifdef TNT_DEBUG + std::cout << "Created Array1D(int n) \n"; +#endif +} + +template +Array1D::Array1D(int n, const T &val) : v_(n), n_(n), data_(v_.begin()) +{ +#ifdef TNT_DEBUG + std::cout << "Created Array1D(int n, const T& val) \n"; +#endif + set_(data_, data_+ n, val); + +} + +template +Array1D::Array1D(int n, T *a) : v_(a), n_(n) , data_(v_.begin()) +{ +#ifdef TNT_DEBUG + std::cout << "Created Array1D(int n, T* a) \n"; +#endif +} + +template +inline Array1D::operator T*() +{ + return &(v_[0]); +} + + +template +inline Array1D::operator const T*() +{ + return &(v_[0]); +} + + + +template +inline T& Array1D::operator[](int i) +{ +#ifdef TNT_BOUNDS_CHECK + assert(i>= 0); + assert(i < n_); +#endif + return data_[i]; +} + +template +inline const T& Array1D::operator[](int i) const +{ +#ifdef TNT_BOUNDS_CHECK + assert(i>= 0); + assert(i < n_); +#endif + return data_[i]; +} + + + + +template +Array1D & Array1D::operator=(const T &a) +{ + set_(data_, data_+n_, a); + return *this; +} + +template +Array1D Array1D::copy() const +{ + Array1D A( n_); + copy_(A.data_, data_, n_); + + return A; +} + + +template +Array1D & Array1D::inject(const Array1D &A) +{ + if (A.n_ == n_) + copy_(data_, A.data_, n_); + + return *this; +} + + + + + +template +Array1D & Array1D::ref(const Array1D &A) +{ + if (this != &A) + { + v_ = A.v_; /* operator= handles the reference counting. */ + n_ = A.n_; + data_ = A.data_; + + } + return *this; +} + +template +Array1D & Array1D::operator=(const Array1D &A) +{ + return ref(A); +} + +template +inline int Array1D::dim1() const { return n_; } + +template +inline int Array1D::dim() const { return n_; } + +template +Array1D::~Array1D() {} + + +/* ............................ exented interface ......................*/ + +template +inline int Array1D::ref_count() const +{ + return v_.ref_count(); +} + +template +inline Array1D Array1D::subarray(int i0, int i1) +{ + if ((i0 > 0) && (i1 < n_) || (i0 <= i1)) + { + Array1D X(*this); /* create a new instance of this array. */ + X.n_ = i1-i0+1; + X.data_ += i0; + + return X; + } + else + { + return Array1D(); + } +} + + +/* private internal functions */ + + +template +void Array1D::set_(T* begin, T* end, const T& a) +{ + for (T* p=begin; p +void Array1D::copy_(T* p, const T* q, int len) const +{ + T *end = p + len; + while (p +#include + +namespace TNT +{ + + +template +std::ostream& operator<<(std::ostream &s, const Array1D &A) +{ + int N=A.dim1(); + +#ifdef TNT_DEBUG + s << "addr: " << (void *) &A[0] << "\n"; +#endif + s << N << "\n"; + for (int j=0; j +std::istream& operator>>(std::istream &s, Array1D &A) +{ + int N; + s >> N; + + Array1D B(N); + for (int i=0; i> B[i]; + A = B; + return s; +} + + + +template +Array1D operator+(const Array1D &A, const Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() != n ) + return Array1D(); + + else + { + Array1D C(n); + + for (int i=0; i +Array1D operator-(const Array1D &A, const Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() != n ) + return Array1D(); + + else + { + Array1D C(n); + + for (int i=0; i +Array1D operator*(const Array1D &A, const Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() != n ) + return Array1D(); + + else + { + Array1D C(n); + + for (int i=0; i +Array1D operator/(const Array1D &A, const Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() != n ) + return Array1D(); + + else + { + Array1D C(n); + + for (int i=0; i +Array1D& operator+=(Array1D &A, const Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() == n) + { + for (int i=0; i +Array1D& operator-=(Array1D &A, const Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() == n) + { + for (int i=0; i +Array1D& operator*=(Array1D &A, const Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() == n) + { + for (int i=0; i +Array1D& operator/=(Array1D &A, const Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() == n) + { + for (int i=0; i +#include +#ifdef TNT_BOUNDS_CHECK +#include +#endif + +#include "tnt_array1d.h" + +namespace TNT +{ + +template +class Array2D +{ + + + private: + + + + Array1D data_; + Array1D v_; + int m_; + int n_; + + public: + + typedef T value_type; + Array2D(); + Array2D(int m, int n); + Array2D(int m, int n, T *a); + Array2D(int m, int n, const T &a); + inline Array2D(const Array2D &A); + inline operator T**(); + inline operator const T**(); + inline Array2D & operator=(const T &a); + inline Array2D & operator=(const Array2D &A); + inline Array2D & ref(const Array2D &A); + Array2D copy() const; + Array2D & inject(const Array2D & A); + inline T* operator[](int i); + inline const T* operator[](int i) const; + inline int dim1() const; + inline int dim2() const; + ~Array2D(); + + /* extended interface (not part of the standard) */ + + + inline int ref_count(); + inline int ref_count_data(); + inline int ref_count_dim1(); + Array2D subarray(int i0, int i1, int j0, int j1); + +}; + + +template +Array2D::Array2D() : data_(), v_(), m_(0), n_(0) {} + +template +Array2D::Array2D(const Array2D &A) : data_(A.data_), v_(A.v_), + m_(A.m_), n_(A.n_) {} + + + + +template +Array2D::Array2D(int m, int n) : data_(m*n), v_(m), m_(m), n_(n) +{ + if (m>0 && n>0) + { + T* p = &(data_[0]); + for (int i=0; i +Array2D::Array2D(int m, int n, const T &val) : data_(m*n), v_(m), + m_(m), n_(n) +{ + if (m>0 && n>0) + { + data_ = val; + T* p = &(data_[0]); + for (int i=0; i +Array2D::Array2D(int m, int n, T *a) : data_(m*n, a), v_(m), m_(m), n_(n) +{ + if (m>0 && n>0) + { + T* p = &(data_[0]); + + for (int i=0; i +inline T* Array2D::operator[](int i) +{ +#ifdef TNT_BOUNDS_CHECK + assert(i >= 0); + assert(i < m_); +#endif + +return v_[i]; + +} + + +template +inline const T* Array2D::operator[](int i) const +{ +#ifdef TNT_BOUNDS_CHECK + assert(i >= 0); + assert(i < m_); +#endif + +return v_[i]; + +} + +template +Array2D & Array2D::operator=(const T &a) +{ + /* non-optimzied, but will work with subarrays in future verions */ + + for (int i=0; i +Array2D Array2D::copy() const +{ + Array2D A(m_, n_); + + for (int i=0; i +Array2D & Array2D::inject(const Array2D &A) +{ + if (A.m_ == m_ && A.n_ == n_) + { + for (int i=0; i +Array2D & Array2D::ref(const Array2D &A) +{ + if (this != &A) + { + v_ = A.v_; + data_ = A.data_; + m_ = A.m_; + n_ = A.n_; + + } + return *this; +} + + + +template +Array2D & Array2D::operator=(const Array2D &A) +{ + return ref(A); +} + +template +inline int Array2D::dim1() const { return m_; } + +template +inline int Array2D::dim2() const { return n_; } + + +template +Array2D::~Array2D() {} + + + + +template +inline Array2D::operator T**() +{ + return &(v_[0]); +} +template +inline Array2D::operator const T**() +{ + return &(v_[0]); +} + +/* ............... extended interface ............... */ +/** + Create a new view to a subarray defined by the boundaries + [i0][i0] and [i1][j1]. The size of the subarray is + (i1-i0) by (j1-j0). If either of these lengths are zero + or negative, the subarray view is null. + +*/ +template +Array2D Array2D::subarray(int i0, int i1, int j0, int j1) +{ + Array2D A; + int m = i1-i0+1; + int n = j1-j0+1; + + /* if either length is zero or negative, this is an invalide + subarray. return a null view. + */ + if (m<1 || n<1) + return A; + + A.data_ = data_; + A.m_ = m; + A.n_ = n; + A.v_ = Array1D(m); + T* p = &(data_[0]) + i0 * n_ + j0; + for (int i=0; i +inline int Array2D::ref_count() +{ + return ref_count_data(); +} + + + +template +inline int Array2D::ref_count_data() +{ + return data_.ref_count(); +} + +template +inline int Array2D::ref_count_dim1() +{ + return v_.ref_count(); +} + + + + +} /* namespace TNT */ + +#endif +/* TNT_ARRAY2D_H */ + diff --git a/lib/include/tnt/tnt_array2d_utils.h b/lib/include/tnt/tnt_array2d_utils.h new file mode 100644 index 0000000..7041ed3 --- /dev/null +++ b/lib/include/tnt/tnt_array2d_utils.h @@ -0,0 +1,287 @@ +/* +* +* Template Numerical Toolkit (TNT) +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + +#ifndef TNT_ARRAY2D_UTILS_H +#define TNT_ARRAY2D_UTILS_H + +#include +#include + +namespace TNT +{ + + +template +std::ostream& operator<<(std::ostream &s, const Array2D &A) +{ + int M=A.dim1(); + int N=A.dim2(); + + s << M << " " << N << "\n"; + + for (int i=0; i +std::istream& operator>>(std::istream &s, Array2D &A) +{ + + int M, N; + + s >> M >> N; + + Array2D B(M,N); + + for (int i=0; i> B[i][j]; + } + + A = B; + return s; +} + + +template +Array2D operator+(const Array2D &A, const Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() != m || B.dim2() != n ) + return Array2D(); + + else + { + Array2D C(m,n); + + for (int i=0; i +Array2D operator-(const Array2D &A, const Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() != m || B.dim2() != n ) + return Array2D(); + + else + { + Array2D C(m,n); + + for (int i=0; i +Array2D operator*(const Array2D &A, const Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() != m || B.dim2() != n ) + return Array2D(); + + else + { + Array2D C(m,n); + + for (int i=0; i +Array2D operator/(const Array2D &A, const Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() != m || B.dim2() != n ) + return Array2D(); + + else + { + Array2D C(m,n); + + for (int i=0; i +Array2D& operator+=(Array2D &A, const Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() == m || B.dim2() == n ) + { + for (int i=0; i +Array2D& operator-=(Array2D &A, const Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() == m || B.dim2() == n ) + { + for (int i=0; i +Array2D& operator*=(Array2D &A, const Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() == m || B.dim2() == n ) + { + for (int i=0; i +Array2D& operator/=(Array2D &A, const Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() == m || B.dim2() == n ) + { + for (int i=0; i +Array2D matmult(const Array2D &A, const Array2D &B) +{ + if (A.dim2() != B.dim1()) + return Array2D(); + + int M = A.dim1(); + int N = A.dim2(); + int K = B.dim2(); + + Array2D C(M,K); + + for (int i=0; i +#include +#ifdef TNT_BOUNDS_CHECK +#include +#endif + +#include "tnt_array1d.h" +#include "tnt_array2d.h" + +namespace TNT +{ + +template +class Array3D +{ + + + private: + Array1D data_; + Array2D v_; + int m_; + int n_; + int g_; + + + public: + + typedef T value_type; + + Array3D(); + Array3D(int m, int n, int g); + Array3D(int m, int n, int g, T val); + Array3D(int m, int n, int g, T *a); + + inline operator T***(); + inline operator const T***(); + inline Array3D(const Array3D &A); + inline Array3D & operator=(const T &a); + inline Array3D & operator=(const Array3D &A); + inline Array3D & ref(const Array3D &A); + Array3D copy() const; + Array3D & inject(const Array3D & A); + + inline T** operator[](int i); + inline const T* const * operator[](int i) const; + inline int dim1() const; + inline int dim2() const; + inline int dim3() const; + ~Array3D(); + + /* extended interface */ + + inline int ref_count(){ return data_.ref_count(); } + Array3D subarray(int i0, int i1, int j0, int j1, + int k0, int k1); +}; + +template +Array3D::Array3D() : data_(), v_(), m_(0), n_(0) {} + +template +Array3D::Array3D(const Array3D &A) : data_(A.data_), + v_(A.v_), m_(A.m_), n_(A.n_), g_(A.g_) +{ +} + + + +template +Array3D::Array3D(int m, int n, int g) : data_(m*n*g), v_(m,n), + m_(m), n_(n), g_(g) +{ + + if (m>0 && n>0 && g>0) + { + T* p = & (data_[0]); + int ng = n_*g_; + + for (int i=0; i +Array3D::Array3D(int m, int n, int g, T val) : data_(m*n*g, val), + v_(m,n), m_(m), n_(n), g_(g) +{ + if (m>0 && n>0 && g>0) + { + + T* p = & (data_[0]); + int ng = n_*g_; + + for (int i=0; i +Array3D::Array3D(int m, int n, int g, T* a) : + data_(m*n*g, a), v_(m,n), m_(m), n_(n), g_(g) +{ + + if (m>0 && n>0 && g>0) + { + T* p = & (data_[0]); + int ng = n_*g_; + + for (int i=0; i +inline T** Array3D::operator[](int i) +{ +#ifdef TNT_BOUNDS_CHECK + assert(i >= 0); + assert(i < m_); +#endif + +return v_[i]; + +} + +template +inline const T* const * Array3D::operator[](int i) const +{ return v_[i]; } + +template +Array3D & Array3D::operator=(const T &a) +{ + for (int i=0; i +Array3D Array3D::copy() const +{ + Array3D A(m_, n_, g_); + for (int i=0; i +Array3D & Array3D::inject(const Array3D &A) +{ + if (A.m_ == m_ && A.n_ == n_ && A.g_ == g_) + + for (int i=0; i +Array3D & Array3D::ref(const Array3D &A) +{ + if (this != &A) + { + m_ = A.m_; + n_ = A.n_; + g_ = A.g_; + v_ = A.v_; + data_ = A.data_; + } + return *this; +} + +template +Array3D & Array3D::operator=(const Array3D &A) +{ + return ref(A); +} + + +template +inline int Array3D::dim1() const { return m_; } + +template +inline int Array3D::dim2() const { return n_; } + +template +inline int Array3D::dim3() const { return g_; } + + + +template +Array3D::~Array3D() {} + +template +inline Array3D::operator T***() +{ + return v_; +} + + +template +inline Array3D::operator const T***() +{ + return v_; +} + +/* extended interface */ +template +Array3D Array3D::subarray(int i0, int i1, int j0, + int j1, int k0, int k1) +{ + + /* check that ranges are valid. */ + if (!( 0 <= i0 && i0 <= i1 && i1 < m_ && + 0 <= j0 && j0 <= j1 && j1 < n_ && + 0 <= k0 && k0 <= k1 && k1 < g_)) + return Array3D(); /* null array */ + + + Array3D A; + A.data_ = data_; + A.m_ = i1-i0+1; + A.n_ = j1-j0+1; + A.g_ = k1-k0+1; + A.v_ = Array2D(A.m_,A.n_); + T* p = &(data_[0]) + i0*n_*g_ + j0*g_ + k0; + + for (int i=0; i +#include + +namespace TNT +{ + + +template +std::ostream& operator<<(std::ostream &s, const Array3D &A) +{ + int M=A.dim1(); + int N=A.dim2(); + int K=A.dim3(); + + s << M << " " << N << " " << K << "\n"; + + for (int i=0; i +std::istream& operator>>(std::istream &s, Array3D &A) +{ + + int M, N, K; + + s >> M >> N >> K; + + Array3D B(M,N,K); + + for (int i=0; i> B[i][j][k]; + + A = B; + return s; +} + + + +template +Array3D operator+(const Array3D &A, const Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) + return Array3D(); + + else + { + Array3D C(m,n,p); + + for (int i=0; i +Array3D operator-(const Array3D &A, const Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) + return Array3D(); + + else + { + Array3D C(m,n,p); + + for (int i=0; i +Array3D operator*(const Array3D &A, const Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) + return Array3D(); + + else + { + Array3D C(m,n,p); + + for (int i=0; i +Array3D operator/(const Array3D &A, const Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) + return Array3D(); + + else + { + Array3D C(m,n,p); + + for (int i=0; i +Array3D& operator+=(Array3D &A, const Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) + { + for (int i=0; i +Array3D& operator-=(Array3D &A, const Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) + { + for (int i=0; i +Array3D& operator*=(Array3D &A, const Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) + { + for (int i=0; i +Array3D& operator/=(Array3D &A, const Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) + { + for (int i=0; i +#include +#include +#include + +namespace TNT +{ + + +template +class Matrix +{ + + + public: + + typedef Subscript size_type; + typedef T value_type; + typedef T element_type; + typedef T* pointer; + typedef T* iterator; + typedef T& reference; + typedef const T* const_iterator; + typedef const T& const_reference; + + Subscript lbound() const { return 1;} + + protected: + Subscript m_; + Subscript n_; + Subscript mn_; // total size + T* v_; + T** row_; + T* vm1_ ; // these point to the same data, but are 1-based + T** rowm1_; + + // internal helper function to create the array + // of row pointers + + void initialize(Subscript M, Subscript N) + { + mn_ = M*N; + m_ = M; + n_ = N; + + v_ = new T[mn_]; + row_ = new T*[M]; + rowm1_ = new T*[M]; + + assert(v_ != NULL); + assert(row_ != NULL); + assert(rowm1_ != NULL); + + T* p = v_; + vm1_ = v_ - 1; + for (Subscript i=0; i &A) + { + initialize(A.m_, A.n_); + copy(A.v_); + } + + Matrix(Subscript M, Subscript N, const T& value = T()) + { + initialize(M,N); + set(value); + } + + Matrix(Subscript M, Subscript N, const T* v) + { + initialize(M,N); + copy(v); + } + + Matrix(Subscript M, Subscript N, const char *s) + { + initialize(M,N); + //std::istrstream ins(s); + std::istringstream ins(s); + + Subscript i, j; + + for (i=0; i> row_[i][j]; + } + + // destructor + // + ~Matrix() + { + destroy(); + } + + + // reallocating + // + Matrix& newsize(Subscript M, Subscript N) + { + if (num_rows() == M && num_cols() == N) + return *this; + + destroy(); + initialize(M,N); + + return *this; + } + + + + + // assignments + // + Matrix& operator=(const Matrix &A) + { + if (v_ == A.v_) + return *this; + + if (m_ == A.m_ && n_ == A.n_) // no need to re-alloc + copy(A.v_); + + else + { + destroy(); + initialize(A.m_, A.n_); + copy(A.v_); + } + + return *this; + } + + Matrix& operator=(const T& scalar) + { + set(scalar); + return *this; + } + + + Subscript dim(Subscript d) const + { +#ifdef TNT_BOUNDS_CHECK + assert( d >= 1); + assert( d <= 2); +#endif + return (d==1) ? m_ : ((d==2) ? n_ : 0); + } + + Subscript num_rows() const { return m_; } + Subscript num_cols() const { return n_; } + + + + + inline T* operator[](Subscript i) + { +#ifdef TNT_BOUNDS_CHECK + assert(0<=i); + assert(i < m_) ; +#endif + return row_[i]; + } + + inline const T* operator[](Subscript i) const + { +#ifdef TNT_BOUNDS_CHECK + assert(0<=i); + assert(i < m_) ; +#endif + return row_[i]; + } + + inline reference operator()(Subscript i) + { +#ifdef TNT_BOUNDS_CHECK + assert(1<=i); + assert(i <= mn_) ; +#endif + return vm1_[i]; + } + + inline const_reference operator()(Subscript i) const + { +#ifdef TNT_BOUNDS_CHECK + assert(1<=i); + assert(i <= mn_) ; +#endif + return vm1_[i]; + } + + + + inline reference operator()(Subscript i, Subscript j) + { +#ifdef TNT_BOUNDS_CHECK + assert(1<=i); + assert(i <= m_) ; + assert(1<=j); + assert(j <= n_); +#endif + return rowm1_[i][j]; + } + + + + inline const_reference operator() (Subscript i, Subscript j) const + { +#ifdef TNT_BOUNDS_CHECK + assert(1<=i); + assert(i <= m_) ; + assert(1<=j); + assert(j <= n_); +#endif + return rowm1_[i][j]; + } + + + + +}; + + +/* *************************** I/O ********************************/ + +template +std::ostream& operator<<(std::ostream &s, const Matrix &A) +{ + Subscript M=A.num_rows(); + Subscript N=A.num_cols(); + + s << M << " " << N << "\n"; + + for (Subscript i=0; i +std::istream& operator>>(std::istream &s, Matrix &A) +{ + + Subscript M, N; + + s >> M >> N; + + if ( !(M == A.num_rows() && N == A.num_cols() )) + { + A.newsize(M,N); + } + + + for (Subscript i=0; i> A[i][j]; + } + + + return s; +} + +// *******************[ basic matrix algorithms ]*************************** + + +template +Matrix operator+(const Matrix &A, + const Matrix &B) +{ + Subscript M = A.num_rows(); + Subscript N = A.num_cols(); + + assert(M==B.num_rows()); + assert(N==B.num_cols()); + + Matrix tmp(M,N); + Subscript i,j; + + for (i=0; i +Matrix operator-(const Matrix &A, + const Matrix &B) +{ + Subscript M = A.num_rows(); + Subscript N = A.num_cols(); + + assert(M==B.num_rows()); + assert(N==B.num_cols()); + + Matrix tmp(M,N); + Subscript i,j; + + for (i=0; i +Matrix mult_element(const Matrix &A, + const Matrix &B) +{ + Subscript M = A.num_rows(); + Subscript N = A.num_cols(); + + assert(M==B.num_rows()); + assert(N==B.num_cols()); + + Matrix tmp(M,N); + Subscript i,j; + + for (i=0; i +Matrix transpose(const Matrix &A) +{ + Subscript M = A.num_rows(); + Subscript N = A.num_cols(); + + Matrix S(N,M); + Subscript i, j; + + for (i=0; i +inline Matrix matmult(const Matrix &A, + const Matrix &B) +{ + +#ifdef TNT_BOUNDS_CHECK + assert(A.num_cols() == B.num_rows()); +#endif + + Subscript M = A.num_rows(); + Subscript N = A.num_cols(); + Subscript K = B.num_cols(); + + Matrix tmp(M,K); + T sum; + + for (Subscript i=0; i +inline Matrix operator*(const Matrix &A, + const Matrix &B) +{ + return matmult(A,B); +} + +template +inline int matmult(Matrix& C, const Matrix &A, + const Matrix &B) +{ + + assert(A.num_cols() == B.num_rows()); + + Subscript M = A.num_rows(); + Subscript N = A.num_cols(); + Subscript K = B.num_cols(); + + C.newsize(M,K); + + T sum; + + const T* row_i; + const T* col_k; + + for (Subscript i=0; i +Vector matmult(const Matrix &A, const Vector &x) +{ + +#ifdef TNT_BOUNDS_CHECK + assert(A.num_cols() == x.dim()); +#endif + + Subscript M = A.num_rows(); + Subscript N = A.num_cols(); + + Vector tmp(M); + T sum; + + for (Subscript i=0; i +inline Vector operator*(const Matrix &A, const Vector &x) +{ + return matmult(A,x); +} + +} // namespace TNT + +#endif +// CMAT_H diff --git a/lib/include/tnt/tnt_fortran_array1d.h b/lib/include/tnt/tnt_fortran_array1d.h new file mode 100644 index 0000000..ad3bba0 --- /dev/null +++ b/lib/include/tnt/tnt_fortran_array1d.h @@ -0,0 +1,267 @@ +/* +* +* Template Numerical Toolkit (TNT) +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + + +#ifndef TNT_FORTRAN_ARRAY1D_H +#define TNT_FORTRAN_ARRAY1D_H + +#include +#include + +#ifdef TNT_BOUNDS_CHECK +#include +#endif + + +#include "tnt_i_refvec.h" + +namespace TNT +{ + +template +class Fortran_Array1D +{ + + private: + + i_refvec v_; + int n_; + T* data_; /* this normally points to v_.begin(), but + * could also point to a portion (subvector) + * of v_. + */ + + void initialize_(int n); + void copy_(T* p, const T* q, int len) const; + void set_(T* begin, T* end, const T& val); + + + public: + + typedef T value_type; + + + Fortran_Array1D(); + explicit Fortran_Array1D(int n); + Fortran_Array1D(int n, const T &a); + Fortran_Array1D(int n, T *a); + inline Fortran_Array1D(const Fortran_Array1D &A); + inline Fortran_Array1D & operator=(const T &a); + inline Fortran_Array1D & operator=(const Fortran_Array1D &A); + inline Fortran_Array1D & ref(const Fortran_Array1D &A); + Fortran_Array1D copy() const; + Fortran_Array1D & inject(const Fortran_Array1D & A); + inline T& operator()(int i); + inline const T& operator()(int i) const; + inline int dim1() const; + inline int dim() const; + ~Fortran_Array1D(); + + + /* ... extended interface ... */ + + inline int ref_count() const; + inline Fortran_Array1D subarray(int i0, int i1); + +}; + + + + +template +Fortran_Array1D::Fortran_Array1D() : v_(), n_(0), data_(0) {} + +template +Fortran_Array1D::Fortran_Array1D(const Fortran_Array1D &A) : v_(A.v_), n_(A.n_), + data_(A.data_) +{ +#ifdef TNT_DEBUG + std::cout << "Created Fortran_Array1D(const Fortran_Array1D &A) \n"; +#endif + +} + + +template +Fortran_Array1D::Fortran_Array1D(int n) : v_(n), n_(n), data_(v_.begin()) +{ +#ifdef TNT_DEBUG + std::cout << "Created Fortran_Array1D(int n) \n"; +#endif +} + +template +Fortran_Array1D::Fortran_Array1D(int n, const T &val) : v_(n), n_(n), data_(v_.begin()) +{ +#ifdef TNT_DEBUG + std::cout << "Created Fortran_Array1D(int n, const T& val) \n"; +#endif + set_(data_, data_+ n, val); + +} + +template +Fortran_Array1D::Fortran_Array1D(int n, T *a) : v_(a), n_(n) , data_(v_.begin()) +{ +#ifdef TNT_DEBUG + std::cout << "Created Fortran_Array1D(int n, T* a) \n"; +#endif +} + +template +inline T& Fortran_Array1D::operator()(int i) +{ +#ifdef TNT_BOUNDS_CHECK + assert(i>= 1); + assert(i <= n_); +#endif + return data_[i-1]; +} + +template +inline const T& Fortran_Array1D::operator()(int i) const +{ +#ifdef TNT_BOUNDS_CHECK + assert(i>= 1); + assert(i <= n_); +#endif + return data_[i-1]; +} + + + + +template +Fortran_Array1D & Fortran_Array1D::operator=(const T &a) +{ + set_(data_, data_+n_, a); + return *this; +} + +template +Fortran_Array1D Fortran_Array1D::copy() const +{ + Fortran_Array1D A( n_); + copy_(A.data_, data_, n_); + + return A; +} + + +template +Fortran_Array1D & Fortran_Array1D::inject(const Fortran_Array1D &A) +{ + if (A.n_ == n_) + copy_(data_, A.data_, n_); + + return *this; +} + + + + + +template +Fortran_Array1D & Fortran_Array1D::ref(const Fortran_Array1D &A) +{ + if (this != &A) + { + v_ = A.v_; /* operator= handles the reference counting. */ + n_ = A.n_; + data_ = A.data_; + + } + return *this; +} + +template +Fortran_Array1D & Fortran_Array1D::operator=(const Fortran_Array1D &A) +{ + return ref(A); +} + +template +inline int Fortran_Array1D::dim1() const { return n_; } + +template +inline int Fortran_Array1D::dim() const { return n_; } + +template +Fortran_Array1D::~Fortran_Array1D() {} + + +/* ............................ exented interface ......................*/ + +template +inline int Fortran_Array1D::ref_count() const +{ + return v_.ref_count(); +} + +template +inline Fortran_Array1D Fortran_Array1D::subarray(int i0, int i1) +{ +#ifdef TNT_DEBUG + std::cout << "entered subarray. \n"; +#endif + if ((i0 > 0) && (i1 < n_) || (i0 <= i1)) + { + Fortran_Array1D X(*this); /* create a new instance of this array. */ + X.n_ = i1-i0+1; + X.data_ += i0; + + return X; + } + else + { +#ifdef TNT_DEBUG + std::cout << "subarray: null return.\n"; +#endif + return Fortran_Array1D(); + } +} + + +/* private internal functions */ + + +template +void Fortran_Array1D::set_(T* begin, T* end, const T& a) +{ + for (T* p=begin; p +void Fortran_Array1D::copy_(T* p, const T* q, int len) const +{ + T *end = p + len; + while (p + +namespace TNT +{ + + +/** + Write an array to a character outstream. Output format is one that can + be read back in via the in-stream operator: one integer + denoting the array dimension (n), followed by n elements, + one per line. + +*/ +template +std::ostream& operator<<(std::ostream &s, const Fortran_Array1D &A) +{ + int N=A.dim1(); + + s << N << "\n"; + for (int j=1; j<=N; j++) + { + s << A(j) << "\n"; + } + s << "\n"; + + return s; +} + +/** + Read an array from a character stream. Input format + is one integer, denoting the dimension (n), followed + by n whitespace-separated elments. Newlines are ignored + +

+ Note: the array being read into references new memory + storage. If the intent is to fill an existing conformant + array, use cin >> B; A.inject(B) ); + instead or read the elements in one-a-time by hand. + + @param s the charater to read from (typically std::in) + @param A the array to read into. +*/ +template +std::istream& operator>>(std::istream &s, Fortran_Array1D &A) +{ + int N; + s >> N; + + Fortran_Array1D B(N); + for (int i=1; i<=N; i++) + s >> B(i); + A = B; + return s; +} + + +template +Fortran_Array1D operator+(const Fortran_Array1D &A, const Fortran_Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() != n ) + return Fortran_Array1D(); + + else + { + Fortran_Array1D C(n); + + for (int i=1; i<=n; i++) + { + C(i) = A(i) + B(i); + } + return C; + } +} + + + +template +Fortran_Array1D operator-(const Fortran_Array1D &A, const Fortran_Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() != n ) + return Fortran_Array1D(); + + else + { + Fortran_Array1D C(n); + + for (int i=1; i<=n; i++) + { + C(i) = A(i) - B(i); + } + return C; + } +} + + +template +Fortran_Array1D operator*(const Fortran_Array1D &A, const Fortran_Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() != n ) + return Fortran_Array1D(); + + else + { + Fortran_Array1D C(n); + + for (int i=1; i<=n; i++) + { + C(i) = A(i) * B(i); + } + return C; + } +} + + +template +Fortran_Array1D operator/(const Fortran_Array1D &A, const Fortran_Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() != n ) + return Fortran_Array1D(); + + else + { + Fortran_Array1D C(n); + + for (int i=1; i<=n; i++) + { + C(i) = A(i) / B(i); + } + return C; + } +} + + + + + + + + + +template +Fortran_Array1D& operator+=(Fortran_Array1D &A, const Fortran_Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() == n) + { + for (int i=1; i<=n; i++) + { + A(i) += B(i); + } + } + return A; +} + + + + +template +Fortran_Array1D& operator-=(Fortran_Array1D &A, const Fortran_Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() == n) + { + for (int i=1; i<=n; i++) + { + A(i) -= B(i); + } + } + return A; +} + + + +template +Fortran_Array1D& operator*=(Fortran_Array1D &A, const Fortran_Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() == n) + { + for (int i=1; i<=n; i++) + { + A(i) *= B(i); + } + } + return A; +} + + + + +template +Fortran_Array1D& operator/=(Fortran_Array1D &A, const Fortran_Array1D &B) +{ + int n = A.dim1(); + + if (B.dim1() == n) + { + for (int i=1; i<=n; i++) + { + A(i) /= B(i); + } + } + return A; +} + + +} // namespace TNT + +#endif diff --git a/lib/include/tnt/tnt_fortran_array2d.h b/lib/include/tnt/tnt_fortran_array2d.h new file mode 100644 index 0000000..f307536 --- /dev/null +++ b/lib/include/tnt/tnt_fortran_array2d.h @@ -0,0 +1,225 @@ +/* +* +* Template Numerical Toolkit (TNT): Two-dimensional Fortran numerical array +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + + +#ifndef TNT_FORTRAN_ARRAY2D_H +#define TNT_FORTRAN_ARRAY2D_H + +#include +#include + +#ifdef TNT_BOUNDS_CHECK +#include +#endif + +#include "tnt_i_refvec.h" + +namespace TNT +{ + +template +class Fortran_Array2D +{ + + + private: + i_refvec v_; + int m_; + int n_; + T* data_; + + + void initialize_(int n); + void copy_(T* p, const T* q, int len); + void set_(T* begin, T* end, const T& val); + + public: + + typedef T value_type; + + Fortran_Array2D(); + Fortran_Array2D(int m, int n); + Fortran_Array2D(int m, int n, T *a); + Fortran_Array2D(int m, int n, const T &a); + inline Fortran_Array2D(const Fortran_Array2D &A); + inline Fortran_Array2D & operator=(const T &a); + inline Fortran_Array2D & operator=(const Fortran_Array2D &A); + inline Fortran_Array2D & ref(const Fortran_Array2D &A); + Fortran_Array2D copy() const; + Fortran_Array2D & inject(const Fortran_Array2D & A); + inline T& operator()(int i, int j); + inline const T& operator()(int i, int j) const ; + inline int dim1() const; + inline int dim2() const; + ~Fortran_Array2D(); + + /* extended interface */ + + inline int ref_count() const; + +}; + +template +Fortran_Array2D::Fortran_Array2D() : v_(), m_(0), n_(0), data_(0) {} + + +template +Fortran_Array2D::Fortran_Array2D(const Fortran_Array2D &A) : v_(A.v_), + m_(A.m_), n_(A.n_), data_(A.data_) {} + + + +template +Fortran_Array2D::Fortran_Array2D(int m, int n) : v_(m*n), m_(m), n_(n), + data_(v_.begin()) {} + +template +Fortran_Array2D::Fortran_Array2D(int m, int n, const T &val) : + v_(m*n), m_(m), n_(n), data_(v_.begin()) +{ + set_(data_, data_+m*n, val); +} + + +template +Fortran_Array2D::Fortran_Array2D(int m, int n, T *a) : v_(a), + m_(m), n_(n), data_(v_.begin()) {} + + + + +template +inline T& Fortran_Array2D::operator()(int i, int j) +{ +#ifdef TNT_BOUNDS_CHECK + assert(i >= 1); + assert(i <= m_); + assert(j >= 1); + assert(j <= n_); +#endif + + return v_[ (j-1)*m_ + (i-1) ]; + +} + +template +inline const T& Fortran_Array2D::operator()(int i, int j) const +{ +#ifdef TNT_BOUNDS_CHECK + assert(i >= 1); + assert(i <= m_); + assert(j >= 1); + assert(j <= n_); +#endif + + return v_[ (j-1)*m_ + (i-1) ]; + +} + + +template +Fortran_Array2D & Fortran_Array2D::operator=(const T &a) +{ + set_(data_, data_+m_*n_, a); + return *this; +} + +template +Fortran_Array2D Fortran_Array2D::copy() const +{ + + Fortran_Array2D B(m_,n_); + + B.inject(*this); + return B; +} + + +template +Fortran_Array2D & Fortran_Array2D::inject(const Fortran_Array2D &A) +{ + if (m_ == A.m_ && n_ == A.n_) + copy_(data_, A.data_, m_*n_); + + return *this; +} + + + +template +Fortran_Array2D & Fortran_Array2D::ref(const Fortran_Array2D &A) +{ + if (this != &A) + { + v_ = A.v_; + m_ = A.m_; + n_ = A.n_; + data_ = A.data_; + } + return *this; +} + +template +Fortran_Array2D & Fortran_Array2D::operator=(const Fortran_Array2D &A) +{ + return ref(A); +} + +template +inline int Fortran_Array2D::dim1() const { return m_; } + +template +inline int Fortran_Array2D::dim2() const { return n_; } + + +template +Fortran_Array2D::~Fortran_Array2D() +{ +} + +template +inline int Fortran_Array2D::ref_count() const { return v_.ref_count(); } + + + + +template +void Fortran_Array2D::set_(T* begin, T* end, const T& a) +{ + for (T* p=begin; p +void Fortran_Array2D::copy_(T* p, const T* q, int len) +{ + T *end = p + len; + while (p + +namespace TNT +{ + + +template +std::ostream& operator<<(std::ostream &s, const Fortran_Array2D &A) +{ + int M=A.dim1(); + int N=A.dim2(); + + s << M << " " << N << "\n"; + + for (int i=1; i<=M; i++) + { + for (int j=1; j<=N; j++) + { + s << A(i,j) << " "; + } + s << "\n"; + } + + + return s; +} + +template +std::istream& operator>>(std::istream &s, Fortran_Array2D &A) +{ + + int M, N; + + s >> M >> N; + + Fortran_Array2D B(M,N); + + for (int i=1; i<=M; i++) + for (int j=1; j<=N; j++) + { + s >> B(i,j); + } + + A = B; + return s; +} + + + + +template +Fortran_Array2D operator+(const Fortran_Array2D &A, const Fortran_Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() != m || B.dim2() != n ) + return Fortran_Array2D(); + + else + { + Fortran_Array2D C(m,n); + + for (int i=1; i<=m; i++) + { + for (int j=1; j<=n; j++) + C(i,j) = A(i,j) + B(i,j); + } + return C; + } +} + +template +Fortran_Array2D operator-(const Fortran_Array2D &A, const Fortran_Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() != m || B.dim2() != n ) + return Fortran_Array2D(); + + else + { + Fortran_Array2D C(m,n); + + for (int i=1; i<=m; i++) + { + for (int j=1; j<=n; j++) + C(i,j) = A(i,j) - B(i,j); + } + return C; + } +} + + +template +Fortran_Array2D operator*(const Fortran_Array2D &A, const Fortran_Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() != m || B.dim2() != n ) + return Fortran_Array2D(); + + else + { + Fortran_Array2D C(m,n); + + for (int i=1; i<=m; i++) + { + for (int j=1; j<=n; j++) + C(i,j) = A(i,j) * B(i,j); + } + return C; + } +} + + +template +Fortran_Array2D operator/(const Fortran_Array2D &A, const Fortran_Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() != m || B.dim2() != n ) + return Fortran_Array2D(); + + else + { + Fortran_Array2D C(m,n); + + for (int i=1; i<=m; i++) + { + for (int j=1; j<=n; j++) + C(i,j) = A(i,j) / B(i,j); + } + return C; + } +} + + + +template +Fortran_Array2D& operator+=(Fortran_Array2D &A, const Fortran_Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() == m || B.dim2() == n ) + { + for (int i=1; i<=m; i++) + { + for (int j=1; j<=n; j++) + A(i,j) += B(i,j); + } + } + return A; +} + +template +Fortran_Array2D& operator-=(Fortran_Array2D &A, const Fortran_Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() == m || B.dim2() == n ) + { + for (int i=1; i<=m; i++) + { + for (int j=1; j<=n; j++) + A(i,j) -= B(i,j); + } + } + return A; +} + +template +Fortran_Array2D& operator*=(Fortran_Array2D &A, const Fortran_Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() == m || B.dim2() == n ) + { + for (int i=1; i<=m; i++) + { + for (int j=1; j<=n; j++) + A(i,j) *= B(i,j); + } + } + return A; +} + +template +Fortran_Array2D& operator/=(Fortran_Array2D &A, const Fortran_Array2D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + + if (B.dim1() == m || B.dim2() == n ) + { + for (int i=1; i<=m; i++) + { + for (int j=1; j<=n; j++) + A(i,j) /= B(i,j); + } + } + return A; +} + +} // namespace TNT + +#endif diff --git a/lib/include/tnt/tnt_fortran_array3d.h b/lib/include/tnt/tnt_fortran_array3d.h new file mode 100644 index 0000000..e51affb --- /dev/null +++ b/lib/include/tnt/tnt_fortran_array3d.h @@ -0,0 +1,223 @@ +/* +* +* Template Numerical Toolkit (TNT): Three-dimensional Fortran numerical array +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + + +#ifndef TNT_FORTRAN_ARRAY3D_H +#define TNT_FORTRAN_ARRAY3D_H + +#include +#include +#ifdef TNT_BOUNDS_CHECK +#include +#endif +#include "tnt_i_refvec.h" + +namespace TNT +{ + +template +class Fortran_Array3D +{ + + + private: + + + i_refvec v_; + int m_; + int n_; + int k_; + T* data_; + + public: + + typedef T value_type; + + Fortran_Array3D(); + Fortran_Array3D(int m, int n, int k); + Fortran_Array3D(int m, int n, int k, T *a); + Fortran_Array3D(int m, int n, int k, const T &a); + inline Fortran_Array3D(const Fortran_Array3D &A); + inline Fortran_Array3D & operator=(const T &a); + inline Fortran_Array3D & operator=(const Fortran_Array3D &A); + inline Fortran_Array3D & ref(const Fortran_Array3D &A); + Fortran_Array3D copy() const; + Fortran_Array3D & inject(const Fortran_Array3D & A); + inline T& operator()(int i, int j, int k); + inline const T& operator()(int i, int j, int k) const ; + inline int dim1() const; + inline int dim2() const; + inline int dim3() const; + inline int ref_count() const; + ~Fortran_Array3D(); + + +}; + +template +Fortran_Array3D::Fortran_Array3D() : v_(), m_(0), n_(0), k_(0), data_(0) {} + + +template +Fortran_Array3D::Fortran_Array3D(const Fortran_Array3D &A) : + v_(A.v_), m_(A.m_), n_(A.n_), k_(A.k_), data_(A.data_) {} + + + +template +Fortran_Array3D::Fortran_Array3D(int m, int n, int k) : + v_(m*n*k), m_(m), n_(n), k_(k), data_(v_.begin()) {} + + + +template +Fortran_Array3D::Fortran_Array3D(int m, int n, int k, const T &val) : + v_(m*n*k), m_(m), n_(n), k_(k), data_(v_.begin()) +{ + for (T* p = data_; p < data_ + m*n*k; p++) + *p = val; +} + +template +Fortran_Array3D::Fortran_Array3D(int m, int n, int k, T *a) : + v_(a), m_(m), n_(n), k_(k), data_(v_.begin()) {} + + + + +template +inline T& Fortran_Array3D::operator()(int i, int j, int k) +{ +#ifdef TNT_BOUNDS_CHECK + assert(i >= 1); + assert(i <= m_); + assert(j >= 1); + assert(j <= n_); + assert(k >= 1); + assert(k <= k_); +#endif + + return data_[(k-1)*m_*n_ + (j-1) * m_ + i-1]; + +} + +template +inline const T& Fortran_Array3D::operator()(int i, int j, int k) const +{ +#ifdef TNT_BOUNDS_CHECK + assert(i >= 1); + assert(i <= m_); + assert(j >= 1); + assert(j <= n_); + assert(k >= 1); + assert(k <= k_); +#endif + + return data_[(k-1)*m_*n_ + (j-1) * m_ + i-1]; +} + + +template +Fortran_Array3D & Fortran_Array3D::operator=(const T &a) +{ + + T *end = data_ + m_*n_*k_; + + for (T *p=data_; p != end; *p++ = a); + + return *this; +} + +template +Fortran_Array3D Fortran_Array3D::copy() const +{ + + Fortran_Array3D B(m_, n_, k_); + B.inject(*this); + return B; + +} + + +template +Fortran_Array3D & Fortran_Array3D::inject(const Fortran_Array3D &A) +{ + + if (m_ == A.m_ && n_ == A.n_ && k_ == A.k_) + { + T *p = data_; + T *end = data_ + m_*n_*k_; + const T* q = A.data_; + for (; p < end; *p++ = *q++); + } + return *this; +} + + + + +template +Fortran_Array3D & Fortran_Array3D::ref(const Fortran_Array3D &A) +{ + + if (this != &A) + { + v_ = A.v_; + m_ = A.m_; + n_ = A.n_; + k_ = A.k_; + data_ = A.data_; + } + return *this; +} + +template +Fortran_Array3D & Fortran_Array3D::operator=(const Fortran_Array3D &A) +{ + return ref(A); +} + +template +inline int Fortran_Array3D::dim1() const { return m_; } + +template +inline int Fortran_Array3D::dim2() const { return n_; } + +template +inline int Fortran_Array3D::dim3() const { return k_; } + + +template +inline int Fortran_Array3D::ref_count() const +{ + return v_.ref_count(); +} + +template +Fortran_Array3D::~Fortran_Array3D() +{ +} + + +} /* namespace TNT */ + +#endif +/* TNT_FORTRAN_ARRAY3D_H */ + diff --git a/lib/include/tnt/tnt_fortran_array3d_utils.h b/lib/include/tnt/tnt_fortran_array3d_utils.h new file mode 100644 index 0000000..a13a275 --- /dev/null +++ b/lib/include/tnt/tnt_fortran_array3d_utils.h @@ -0,0 +1,249 @@ +/* +* +* Template Numerical Toolkit (TNT) +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + +#ifndef TNT_FORTRAN_ARRAY3D_UTILS_H +#define TNT_FORTRAN_ARRAY3D_UTILS_H + +#include +#include + +namespace TNT +{ + + +template +std::ostream& operator<<(std::ostream &s, const Fortran_Array3D &A) +{ + int M=A.dim1(); + int N=A.dim2(); + int K=A.dim3(); + + s << M << " " << N << " " << K << "\n"; + + for (int i=1; i<=M; i++) + { + for (int j=1; j<=N; j++) + { + for (int k=1; k<=K; k++) + s << A(i,j,k) << " "; + s << "\n"; + } + s << "\n"; + } + + + return s; +} + +template +std::istream& operator>>(std::istream &s, Fortran_Array3D &A) +{ + + int M, N, K; + + s >> M >> N >> K; + + Fortran_Array3D B(M,N,K); + + for (int i=1; i<=M; i++) + for (int j=1; j<=N; j++) + for (int k=1; k<=K; k++) + s >> B(i,j,k); + + A = B; + return s; +} + + +template +Fortran_Array3D operator+(const Fortran_Array3D &A, const Fortran_Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) + return Fortran_Array3D(); + + else + { + Fortran_Array3D C(m,n,p); + + for (int i=1; i<=m; i++) + for (int j=1; j<=n; j++) + for (int k=1; k<=p; k++) + C(i,j,k) = A(i,j,k)+ B(i,j,k); + + return C; + } +} + + +template +Fortran_Array3D operator-(const Fortran_Array3D &A, const Fortran_Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) + return Fortran_Array3D(); + + else + { + Fortran_Array3D C(m,n,p); + + for (int i=1; i<=m; i++) + for (int j=1; j<=n; j++) + for (int k=1; k<=p; k++) + C(i,j,k) = A(i,j,k)- B(i,j,k); + + return C; + } +} + + +template +Fortran_Array3D operator*(const Fortran_Array3D &A, const Fortran_Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) + return Fortran_Array3D(); + + else + { + Fortran_Array3D C(m,n,p); + + for (int i=1; i<=m; i++) + for (int j=1; j<=n; j++) + for (int k=1; k<=p; k++) + C(i,j,k) = A(i,j,k)* B(i,j,k); + + return C; + } +} + + +template +Fortran_Array3D operator/(const Fortran_Array3D &A, const Fortran_Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() != m || B.dim2() != n || B.dim3() != p ) + return Fortran_Array3D(); + + else + { + Fortran_Array3D C(m,n,p); + + for (int i=1; i<=m; i++) + for (int j=1; j<=n; j++) + for (int k=1; k<=p; k++) + C(i,j,k) = A(i,j,k)/ B(i,j,k); + + return C; + } +} + + +template +Fortran_Array3D& operator+=(Fortran_Array3D &A, const Fortran_Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) + { + for (int i=1; i<=m; i++) + for (int j=1; j<=n; j++) + for (int k=1; k<=p; k++) + A(i,j,k) += B(i,j,k); + } + + return A; +} + + +template +Fortran_Array3D& operator-=(Fortran_Array3D &A, const Fortran_Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) + { + for (int i=1; i<=m; i++) + for (int j=1; j<=n; j++) + for (int k=1; k<=p; k++) + A(i,j,k) -= B(i,j,k); + } + + return A; +} + + +template +Fortran_Array3D& operator*=(Fortran_Array3D &A, const Fortran_Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) + { + for (int i=1; i<=m; i++) + for (int j=1; j<=n; j++) + for (int k=1; k<=p; k++) + A(i,j,k) *= B(i,j,k); + } + + return A; +} + + +template +Fortran_Array3D& operator/=(Fortran_Array3D &A, const Fortran_Array3D &B) +{ + int m = A.dim1(); + int n = A.dim2(); + int p = A.dim3(); + + if (B.dim1() == m && B.dim2() == n && B.dim3() == p ) + { + for (int i=1; i<=m; i++) + for (int j=1; j<=n; j++) + for (int k=1; k<=p; k++) + A(i,j,k) /= B(i,j,k); + } + + return A; +} + + +} // namespace TNT + +#endif diff --git a/lib/include/tnt/tnt_i_refvec.h b/lib/include/tnt/tnt_i_refvec.h new file mode 100644 index 0000000..5a67eb5 --- /dev/null +++ b/lib/include/tnt/tnt_i_refvec.h @@ -0,0 +1,243 @@ +/* +* +* Template Numerical Toolkit (TNT) +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + + +#ifndef TNT_I_REFVEC_H +#define TNT_I_REFVEC_H + +#include +#include + +#ifdef TNT_BOUNDS_CHECK +#include +#endif + +#ifndef NULL +#define NULL 0 +#endif + +namespace TNT +{ +/* + Internal representation of ref-counted array. The TNT + arrays all use this building block. + +

+ If an array block is created by TNT, then every time + an assignment is made, the left-hand-side reference + is decreased by one, and the right-hand-side refernce + count is increased by one. If the array block was + external to TNT, the refernce count is a NULL pointer + regardless of how many references are made, since the + memory is not freed by TNT. + + + +*/ +template +class i_refvec +{ + + + private: + T* data_; + int *ref_count_; + + + public: + + i_refvec(); + explicit i_refvec(int n); + inline i_refvec(T* data); + inline i_refvec(const i_refvec &v); + inline T* begin(); + inline const T* begin() const; + inline T& operator[](int i); + inline const T& operator[](int i) const; + inline i_refvec & operator=(const i_refvec &V); + void copy_(T* p, const T* q, const T* e); + void set_(T* p, const T* b, const T* e); + inline int ref_count() const; + inline int is_null() const; + inline void destroy(); + ~i_refvec(); + +}; + +template +void i_refvec::copy_(T* p, const T* q, const T* e) +{ + for (T* t=p; q +i_refvec::i_refvec() : data_(NULL), ref_count_(NULL) {} + +/** + In case n is 0 or negative, it does NOT call new. +*/ +template +i_refvec::i_refvec(int n) : data_(NULL), ref_count_(NULL) +{ + if (n >= 1) + { +#ifdef TNT_DEBUG + std::cout << "new data storage.\n"; +#endif + data_ = new T[n]; + ref_count_ = new int; + *ref_count_ = 1; + } +} + +template +inline i_refvec::i_refvec(const i_refvec &V): data_(V.data_), + ref_count_(V.ref_count_) +{ + if (V.ref_count_ != NULL) + (*(V.ref_count_))++; +} + + +template +i_refvec::i_refvec(T* data) : data_(data), ref_count_(NULL) {} + +template +inline T* i_refvec::begin() +{ + return data_; +} + +template +inline const T& i_refvec::operator[](int i) const +{ + return data_[i]; +} + +template +inline T& i_refvec::operator[](int i) +{ + return data_[i]; +} + + +template +inline const T* i_refvec::begin() const +{ + return data_; +} + + + +template +i_refvec & i_refvec::operator=(const i_refvec &V) +{ + if (this == &V) + return *this; + + + if (ref_count_ != NULL) + { + (*ref_count_) --; + if ((*ref_count_) == 0) + destroy(); + } + + data_ = V.data_; + ref_count_ = V.ref_count_; + + if (V.ref_count_ != NULL) + (*(V.ref_count_))++; + + return *this; +} + +template +void i_refvec::destroy() +{ + if (ref_count_ != NULL) + { +#ifdef TNT_DEBUG + std::cout << "destorying data... \n"; +#endif + delete ref_count_; + +#ifdef TNT_DEBUG + std::cout << "deleted ref_count_ ...\n"; +#endif + if (data_ != NULL) + delete []data_; +#ifdef TNT_DEBUG + std::cout << "deleted data_[] ...\n"; +#endif + data_ = NULL; + } +} + +/* +* return 1 is vector is empty, 0 otherwise +* +* if is_null() is false and ref_count() is 0, then +* +*/ +template +int i_refvec::is_null() const +{ + return (data_ == NULL ? 1 : 0); +} + +/* +* returns -1 if data is external, +* returns 0 if a is NULL array, +* otherwise returns the positive number of vectors sharing +* this data space. +*/ +template +int i_refvec::ref_count() const +{ + if (data_ == NULL) + return 0; + else + return (ref_count_ != NULL ? *ref_count_ : -1) ; +} + +template +i_refvec::~i_refvec() +{ + if (ref_count_ != NULL) + { + (*ref_count_)--; + + if (*ref_count_ == 0) + destroy(); + } +} + + +} /* namespace TNT */ + + + + + +#endif +/* TNT_I_REFVEC_H */ + diff --git a/lib/include/tnt/tnt_math_utils.h b/lib/include/tnt/tnt_math_utils.h new file mode 100644 index 0000000..f9c1c91 --- /dev/null +++ b/lib/include/tnt/tnt_math_utils.h @@ -0,0 +1,34 @@ +#ifndef MATH_UTILS_H +#define MATH_UTILS_H + +/* needed for fabs, sqrt() below */ +#include + + + +namespace TNT +{ +/** + @returns hypotenuse of real (non-complex) scalars a and b by + avoiding underflow/overflow + using (a * sqrt( 1 + (b/a) * (b/a))), rather than + sqrt(a*a + b*b). +*/ +template +Real hypot(const Real &a, const Real &b) +{ + + if (a== 0) + return abs(b); + else + { + Real c = b/a; + return fabs(a) * sqrt(1 + c*c); + } +} +} /* TNT namespace */ + + + +#endif +/* MATH_UTILS_H */ diff --git a/lib/include/tnt/tnt_sparse_matrix_csr.h b/lib/include/tnt/tnt_sparse_matrix_csr.h new file mode 100644 index 0000000..0d4fde1 --- /dev/null +++ b/lib/include/tnt/tnt_sparse_matrix_csr.h @@ -0,0 +1,103 @@ +/* +* +* Template Numerical Toolkit (TNT) +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + +#ifndef TNT_SPARSE_MATRIX_CSR_H +#define TNT_SPARSE_MATRIX_CSR_H + +#include "tnt_array1d.h" + +namespace TNT +{ + + +/** + Read-only view of a sparse matrix in compressed-row storage + format. Neither array elements (nonzeros) nor sparsity + structure can be modified. If modifications are required, + create a new view. + +

+ Index values begin at 0. + +

+ Storage requirements: An (m x n) matrix with + nz nonzeros requires no more than ((T+I)*nz + M*I) + bytes, where T is the size of data elements and + I is the size of integers. + + +*/ +template +class Sparse_Matrix_CompRow { + +private: + Array1D val_; // data values (nz_ elements) + Array1D rowptr_; // row_ptr (dim_[0]+1 elements) + Array1D colind_; // col_ind (nz_ elements) + + int dim1_; // number of rows + int dim2_; // number of cols + +public: + + Sparse_Matrix_CompRow(const Sparse_Matrix_CompRow &S); + Sparse_Matrix_CompRow(int M, int N, int nz, const T *val, + const int *r, const int *c); + + + + inline const T& val(int i) const { return val_[i]; } + inline const int& row_ptr(int i) const { return rowptr_[i]; } + inline const int& col_ind(int i) const { return colind_[i];} + + inline int dim1() const {return dim1_;} + inline int dim2() const {return dim2_;} + int NumNonzeros() const {return val_.dim1();} + + + Sparse_Matrix_CompRow& operator=( + const Sparse_Matrix_CompRow &R); + + + +}; + +/** + Construct a read-only view of existing sparse matrix in + compressed-row storage format. + + @param M the number of rows of sparse matrix + @param N the number of columns of sparse matrix + @param nz the number of nonzeros + @param val a contiguous list of nonzero values + @param r row-pointers: r[i] denotes the begining position of row i + (i.e. the ith row begins at val[row[i]]). + @param c column-indices: c[i] denotes the column location of val[i] +*/ +template +Sparse_Matrix_CompRow::Sparse_Matrix_CompRow(int M, int N, int nz, + const T *val, const int *r, const int *c) : val_(nz,val), + rowptr_(M, r), colind_(nz, c), dim1_(M), dim2_(N) {} + + +} +// namespace TNT + +#endif diff --git a/lib/include/tnt/tnt_stopwatch.h b/lib/include/tnt/tnt_stopwatch.h new file mode 100644 index 0000000..8dc5d23 --- /dev/null +++ b/lib/include/tnt/tnt_stopwatch.h @@ -0,0 +1,95 @@ +/* +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + + +#ifndef STOPWATCH_H +#define STOPWATCH_H + +// for clock() and CLOCKS_PER_SEC +#include + + +namespace TNT +{ + +inline static double seconds(void) +{ + const double secs_per_tick = 1.0 / CLOCKS_PER_SEC; + return ( (double) clock() ) * secs_per_tick; +} + +class Stopwatch { + private: + int running_; + double start_time_; + double total_; + + public: + inline Stopwatch(); + inline void start(); + inline double stop(); + inline double read(); + inline void resume(); + inline int running(); +}; + +inline Stopwatch::Stopwatch() : running_(0), start_time_(0.0), total_(0.0) {} + +void Stopwatch::start() +{ + running_ = 1; + total_ = 0.0; + start_time_ = seconds(); +} + +double Stopwatch::stop() +{ + if (running_) + { + total_ += (seconds() - start_time_); + running_ = 0; + } + return total_; +} + +inline void Stopwatch::resume() +{ + if (!running_) + { + start_time_ = seconds(); + running_ = 1; + } +} + + +inline double Stopwatch::read() +{ + if (running_) + { + stop(); + resume(); + } + return total_; +} + + +} /* TNT namespace */ +#endif + + + diff --git a/lib/include/tnt/tnt_subscript.h b/lib/include/tnt/tnt_subscript.h new file mode 100644 index 0000000..d8fe120 --- /dev/null +++ b/lib/include/tnt/tnt_subscript.h @@ -0,0 +1,54 @@ +/* +* +* Template Numerical Toolkit (TNT) +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + +#ifndef TNT_SUBSCRPT_H +#define TNT_SUBSCRPT_H + + +//--------------------------------------------------------------------- +// This definition describes the default TNT data type used for +// indexing into TNT matrices and vectors. The data type should +// be wide enough to index into large arrays. It defaults to an +// "int", but can be overriden at compile time redefining TNT_SUBSCRIPT_TYPE, +// e.g. +// +// c++ -DTNT_SUBSCRIPT_TYPE='unsigned int' ... +// +//--------------------------------------------------------------------- +// + +#ifndef TNT_SUBSCRIPT_TYPE +#define TNT_SUBSCRIPT_TYPE int +#endif + +namespace TNT +{ + typedef TNT_SUBSCRIPT_TYPE Subscript; +} /* namespace TNT */ + + +// () indexing in TNT means 1-offset, i.e. x(1) and A(1,1) are the +// first elements. This offset is left as a macro for future +// purposes, but should not be changed in the current release. +// +// +#define TNT_BASE_OFFSET (1) + +#endif diff --git a/lib/include/tnt/tnt_vec.h b/lib/include/tnt/tnt_vec.h new file mode 100644 index 0000000..3455d79 --- /dev/null +++ b/lib/include/tnt/tnt_vec.h @@ -0,0 +1,404 @@ +/* +* +* Template Numerical Toolkit (TNT) +* +* Mathematical and Computational Sciences Division +* National Institute of Technology, +* Gaithersburg, MD USA +* +* +* This software was developed at the National Institute of Standards and +* Technology (NIST) by employees of the Federal Government in the course +* of their official duties. Pursuant to title 17 Section 105 of the +* United States Code, this software is not subject to copyright protection +* and is in the public domain. NIST assumes no responsibility whatsoever for +* its use by other parties, and makes no guarantees, expressed or implied, +* about its quality, reliability, or any other characteristic. +* +*/ + + + +#ifndef TNT_VEC_H +#define TNT_VEC_H + +#include "tnt_subscript.h" +#include +#include +#include +#include + +namespace TNT +{ + +/** + [Deprecatred] Value-based vector class from pre-1.0 + TNT version. Kept here for backward compatiblity, but should + use the newer TNT::Array1D classes instead. + +*/ + +template +class Vector +{ + + + public: + + typedef Subscript size_type; + typedef T value_type; + typedef T element_type; + typedef T* pointer; + typedef T* iterator; + typedef T& reference; + typedef const T* const_iterator; + typedef const T& const_reference; + + Subscript lbound() const { return 1;} + + protected: + T* v_; + T* vm1_; // pointer adjustment for optimzied 1-offset indexing + Subscript n_; + + // internal helper function to create the array + // of row pointers + + void initialize(Subscript N) + { + // adjust pointers so that they are 1-offset: + // v_[] is the internal contiguous array, it is still 0-offset + // + assert(v_ == NULL); + v_ = new T[N]; + assert(v_ != NULL); + vm1_ = v_-1; + n_ = N; + } + + void copy(const T* v) + { + Subscript N = n_; + Subscript i; + +#ifdef TNT_UNROLL_LOOPS + Subscript Nmod4 = N & 3; + Subscript N4 = N - Nmod4; + + for (i=0; i &A) : v_(0), vm1_(0), n_(0) + { + initialize(A.n_); + copy(A.v_); + } + + Vector(Subscript N, const T& value = T()) : v_(0), vm1_(0), n_(0) + { + initialize(N); + set(value); + } + + Vector(Subscript N, const T* v) : v_(0), vm1_(0), n_(0) + { + initialize(N); + copy(v); + } + + Vector(Subscript N, char *s) : v_(0), vm1_(0), n_(0) + { + initialize(N); + std::istringstream ins(s); + + Subscript i; + + for (i=0; i> v_[i]; + } + + + // methods + // + Vector& newsize(Subscript N) + { + if (n_ == N) return *this; + + destroy(); + initialize(N); + + return *this; + } + + + // assignments + // + Vector& operator=(const Vector &A) + { + if (v_ == A.v_) + return *this; + + if (n_ == A.n_) // no need to re-alloc + copy(A.v_); + + else + { + destroy(); + initialize(A.n_); + copy(A.v_); + } + + return *this; + } + + Vector& operator=(const T& scalar) + { + set(scalar); + return *this; + } + + inline Subscript dim() const + { + return n_; + } + + inline Subscript size() const + { + return n_; + } + + + inline reference operator()(Subscript i) + { +#ifdef TNT_BOUNDS_CHECK + assert(1<=i); + assert(i <= n_) ; +#endif + return vm1_[i]; + } + + inline const_reference operator() (Subscript i) const + { +#ifdef TNT_BOUNDS_CHECK + assert(1<=i); + assert(i <= n_) ; +#endif + return vm1_[i]; + } + + inline reference operator[](Subscript i) + { +#ifdef TNT_BOUNDS_CHECK + assert(0<=i); + assert(i < n_) ; +#endif + return v_[i]; + } + + inline const_reference operator[](Subscript i) const + { +#ifdef TNT_BOUNDS_CHECK + assert(0<=i); + + + + + + + assert(i < n_) ; +#endif + return v_[i]; + } + + + +}; + + +/* *************************** I/O ********************************/ + +template +std::ostream& operator<<(std::ostream &s, const Vector &A) +{ + Subscript N=A.dim(); + + s << N << "\n"; + + for (Subscript i=0; i +std::istream & operator>>(std::istream &s, Vector &A) +{ + + Subscript N; + + s >> N; + + if ( !(N == A.size() )) + { + A.newsize(N); + } + + + for (Subscript i=0; i> A[i]; + + + return s; +} + +// *******************[ basic matrix algorithms ]*************************** + + +template +Vector operator+(const Vector &A, + const Vector &B) +{ + Subscript N = A.dim(); + + assert(N==B.dim()); + + Vector tmp(N); + Subscript i; + + for (i=0; i +Vector operator-(const Vector &A, + const Vector &B) +{ + Subscript N = A.dim(); + + assert(N==B.dim()); + + Vector tmp(N); + Subscript i; + + for (i=0; i +Vector operator*(const Vector &A, + const Vector &B) +{ + Subscript N = A.dim(); + + assert(N==B.dim()); + + Vector tmp(N); + Subscript i; + + for (i=0; i +T dot_prod(const Vector &A, const Vector &B) +{ + Subscript N = A.dim(); + assert(N == B.dim()); + + Subscript i; + T sum = 0; + + for (i=0; i