diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..51400c0 --- /dev/null +++ b/Makefile @@ -0,0 +1,52 @@ +# application name +APPNAME := gtksolver +# compiler +CC := gcc +CCLD := $(CC) +# output/object/include/source directories +BINDIR := bin +OBJDIR := obj +INCLUDE := include +SRCDIR := src +# compiler and linker flags +CFLAGS := -Wall -Wextra -pedantic -finline-functions -std=c11 -Wshadow +CFLAGS += -I$(INCLUDE) +CFLAGS += `pkg-config --cflags --libs gtk+-2.0` +ifeq ($(debug),-DDEBUG) +CFLAGS += -g +else +CFLAGS += -Ofast +endif +LDFLAGS := `pkg-config --libs gtk+-2.0` +# libraries +LIBS := +# source/include/object variables +SOURCES := $(wildcard $(SRCDIR)/*.c) +INCLUDES := $(wildcard $(INCLUDE)/*.h) +OBJECTS := $(SOURCES:$(SRCDIR)/%.c=$(OBJDIR)/%.o) + +ifeq ($(os),windows) +APPNAME := $(APPNAME).exe +LIBS := -Wl,-subsystem,windows +endif + +# debug printing variable contents +# $(info $$SOURCES is [${SOURCES}]) +# $(info $$INCLUDES is [${INCLUDES}]) +# $(info $$OBJECTS is [${OBJECTS}]) + +# $(APPNAME): $(OBJECTS) +all: $(OBJECTS) + @mkdir -p $(@D)/bin + $(CCLD) -o $(BINDIR)/$(APPNAME) $(OBJECTS) $(CFLAGS) $(LDFLAGS) $(LIBS) +# strip only if -DDEBUG not set +ifneq ($(debug),-DDEBUG) + strip -s $(BINDIR)/$(APPNAME) +endif + +$(OBJECTS): $(OBJDIR)/%.o : $(SRCDIR)/%.c + @mkdir -p $(OBJDIR) + $(CC) $(CFLAGS) -c -o $@ $< + +clean: + rm -rf $(BINDIR) $(OBJDIR) diff --git a/README.md b/README.md new file mode 100644 index 0000000..76e7c62 --- /dev/null +++ b/README.md @@ -0,0 +1,85 @@ +# gtksolver +**GtkSolver** - Linear System of Equation Solver written in C and Gtk+2. + +GtkSolver is a simple but fast and capable linear system solver that is simply a GtkTextView wrapped around a command-line solver to allow simple copy/paste of the coefficent matrix (including the constant vector as the last column) into the textview editor for solving. Basically it is an editor with a `[Solve...]` button. A convenience that saves creating the file with the coefficent matrix and constant vector before calling the solver from the command line. The underlying parser and solver is written entirely in C. (the `mtrx_t.[ch]` source files contain code for handling direct file input as well) + +The linear systems solver uses Gauss-Jordan Elimination with full pivoting to solve a system of equations contianing any number of unknowns up to the physical memory limits of your computer. See [Gaussian Elimination](https://en.wikipedia.org/wiki/Gaussian_elimination) + +### Solver Use + +The interface is straight forward. The program lauches with a short help message and example of the input format expected in the textview itself. Simply paste your coefficent matrix with the constant vector as the final column over the help message and click `[Solve...]`. Or if you need to middle-mouse paste from your select-buffer, the `[Clear]` button allows you to clear the textview and paste without disturbing the contents of your select buffer. + +The contents of the textview (i.e. the GtkTextBuffer) is read and passed to the solver. The input format is flexible, but must be an `[N x N+1]` matrix (representing `N` equations and `N` unknowns PLUS the constant vector as the last column). Delimiters between the values are ignore. The all values are processed as type `double`. (as defined in `mtrx_t.h`) + +The contents of the textview before pressing `[Solve...]` could equally be: + + 3.0 2.0 -4.0 3.0 + 2.0 3.0 3.0 15.0 + 5.0 -3.0 1.0 14.0 + +or + + + linear system of equations: + + 3 2 -4 | 3 + 2 3 3 | 15 + 5 -3 1 | 14 + +or + + 3,2,-4,3 + 2,3,3,15 + 5,-3,1,14 + +(all leaning characters that are not `.+-[0-9]` are discarded, so the line and leading whitespace in the second example above `"linear system of equations:"` is ignored) + +Clicking `[Clear]` clears the GtkTextBuffer, and clicking `[Help]` clears the buffer and redisplays the initial help message. + +### Output + +Taking the contents of the second example above as the contents of the textview and clicking `[Solve...]` results in: + + linear system of equations: + + 3 2 -4 | 3 + 2 3 3 | 15 + 5 -3 1 | 14 + + + Solution Vector: + + x[ 0] : 3.0000000 + x[ 1] : 1.0000000 + x[ 2] : 2.0000000 + +(where the formatted solution vector is simply written back to the same GtkTextBuffer and displayed in the textview) + +### Compiling + +For Linux, all that is needed is `gcc/make/pkg-config` and `Gtk+2`. (note some distributions package the headers and development files in separate packages, for instance `Gtk+2-dev`). You may want to create an out-of-source directory for building to prevent cluttering your sources with the object and executable files. Simply create a separate directory (e.g. `gtksolver.build`) and then symlink the `Makefile`, `src` and `include` directories. All that is needed then is to change to the build directory and type: + + $ make + +For building on Windows, see the notes on obtaining the precompiled Gtk libraries and header files in the [GtkWrite Readme.md](https://github.com/drankinatty/gtkwrite). + +Building with **MinGW** on windows: + + $ make os=windows + +(the `os=windows` simply appends `.exe` to the executable name so that `strip` finds it correctly.) + + +Building with **TDM-MinGW** on windows. + +The only addition for using TDM-MinGW is you must set the `CC=mingw32-gcc` make-variable and invoke make with `mingw32-make` given the naming difference, e.g. + + $ mingw32-make CC=mingw32-gcc os=windows + +### Development Status + +As mentioned at the beginning, this project is basically a quick GtkTextView wrapper around a linear system solver. The Gtk+2 interface and fowarding of the textview contents to the solver has minimal validations. This basically grew out of finding a convenient way to help students with physics and engineering problems. The underlying solver and parsing code is much more robust. + +### License/Copyright + +This code is licensed under the GPLv2 open-source license contained in `gpl-2.0.txt` and copyright David C. Rankin, 2019. diff --git a/gpl-2.0.txt b/gpl-2.0.txt new file mode 100644 index 0000000..d159169 --- /dev/null +++ b/gpl-2.0.txt @@ -0,0 +1,339 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Lesser General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. diff --git a/include/memrealloc.h b/include/memrealloc.h new file mode 100644 index 0000000..38a37b8 --- /dev/null +++ b/include/memrealloc.h @@ -0,0 +1,29 @@ +#ifndef _memrealloc_h_ +#define _memrealloc_h_ 1 + +#include +#include +#include + +/** realloc 'ptr' of 'nelem' of 'psz' from 'oldsz' to 'newsz' of 'psz'. + * returns pointer to reallocated block of memory with new + * memory initialized to 0/NULL. return must be assigned to + * original pointer in caller. + */ +void *xrealloc_fixed (void *ptr, size_t psz, size_t oldsz, size_t newsz); + +/** realloc 'ptr' of 'nelem' of 'psz' to 'nelem + inc' of 'psz'. + * returns pointer to reallocated block of memory with new + * memory initialized to 0/NULL. return must be assigned to + * original pointer in caller. + */ +void *xrealloc_inc (void *ptr, size_t psz, size_t *nelem, size_t inc); + +/** realloc 'ptr' of 'nelem' of 'psz' to 'nelem * 2' of 'psz'. + * returns pointer to reallocated block of memory with new + * memory initialized to 0/NULL. return must be assigned to + * original pointer in caller. + */ +void *xrealloc_x2 (void *ptr, size_t psz, size_t *nelem); + +#endif diff --git a/include/mtrx_t.h b/include/mtrx_t.h new file mode 100644 index 0000000..276528f --- /dev/null +++ b/include/mtrx_t.h @@ -0,0 +1,125 @@ +#ifndef __mtrx_t_h__ +#define __mtrx_t_h__ 1 + +#include +#include +#include +#include +#include +#include + +typedef double T; /* must set T typedef for type */ + +typedef struct { /* matrix struct */ + size_t rows, cols, rowmax, colmax; + T **mtrx; +} mtrx_t; + +typedef struct { /* vector struct (not fully implemented) */ + size_t nelem, nmax; + T *vect; +} vect_t; + +#define ROWSZ 8 /* initial allocation size */ +#define COLSZ 8 + +#define MAXC 1024 /* max chars per-row to read */ + +#define T_FP_MAX 1.0e-20 /* floating-point values considered zero */ +#define T_FP_MIN -1.0e-20 + +enum { ROWOP, COLOP }; /* row operation/column operation consts */ + +/* __func__ is not defined in pre C99 so provide fallback */ +#if __STDC_VERSION__ < 199901L +# if __GNUC__ >= 2 +# define __func__ __FUNCTION__ +# else +# define __func__ "" +# endif +#endif + +/** allocate vector of n elements */ +T *vect_calloc (const size_t n); +/** allocate matrix of (m x n) */ +T **mtrx_calloc (const size_t m, const size_t n); +/** allocate/initialize struct T and (m) pointers, set rowmax/colmax */ +mtrx_t *mtrx_create_ptrs (const size_t m); +/** allocate/initialize struct T and (m x n) matrix, set rowmax/colmax */ +mtrx_t *mtrx_create_fixed (const size_t m, const size_t n); +/** realloc matrix pointers to twice current */ +T **mtrx_realloc_ptrs (mtrx_t *m); +/** realloc matrix_fixed, which = 0 realloc ptrs, 1 column storage. + * custom realloc used as no need to zero new pointers (they will be + * allocated) and row realloc simply saves the function call overhead + * for each row realloc. + */ +T **mtrx_realloc_fixed (mtrx_t *m, const int which); +/** read (m x n) matrix from file stream into dynamically sized T *mtrx. + * pointers are reallocated x2 as needed, row storage is reallocated + * x2 as required to store all values in 1st row of data, then row storage + * is reallocate to shrink row allocation to exact number of columns and + * T->colmax is set. following read of all data, pointer storage is reduced + * to T->nrows and T->rowmax is set. + */ +mtrx_t *mtrx_read_alloc (FILE *fp); +mtrx_t *mtrx_read_alloc_buf (char *textviewbuf); +/** matr_read_fixed, alloc/read (m x n) matrix from file stream. + * reallocate as required, but should not be required if correct initial + * size given. matrix is resized to exact number of rows and pointers + * before return. + */ +mtrx_t *mtrx_read_fixed (FILE *fp, const size_t rows, const size_t cols); +/** output matrix to stdout with/pad */ +void mtrx_prn (const mtrx_t *m, int width); +/** print a (m x m+1) system of linear equations with separtor */ +void mtrx_sys_prn (const mtrx_t *m, int width); +/** print a (m x m) matrix from (m x n) with n > m */ +void mtrx_prn_sq (const mtrx_t *m, int width); +/** copy a struct matrix */ +mtrx_t *mtrx_copy (const mtrx_t *ma); +/** free a (m x n) matrix and struct for m->nrows pointers */ +void mtrx_free (mtrx_t *m); + +/** add two struct matricies */ +mtrx_t *mtrx_add (const mtrx_t *ma, const mtrx_t *mb); +/** subtract two struct matricies */ +mtrx_t *mtrx_sub (const mtrx_t *ma, const mtrx_t *mb); +/** multiply two struct matricies */ +mtrx_t *mtrx_mult (const mtrx_t *ma, const mtrx_t *mb); +/** transpose matrix */ +mtrx_t *mtrx_trans (const mtrx_t *m); + +/** solve system of equations */ +T *mtrx_solv (const mtrx_t *m, const T *v); +/** mtrx_solv_cmb - solves system of eq where m contains solution vect. + * solves system where m contains the solution vector as the last column + * in the form {m} = {m'}[v] where {m'} is the coefficient matrix and [v] + * the solution vector. it is a wrapper for mtrx_solv() above. returns + * a vector containing the unique solution or NULL if {m'} is singular + * or the solutions are infinite or trivial. + */ +T *mtrx_solv_cmb (mtrx_t *m); +/** Guass-Jordan elimination with full pivoting. + * 'a' is coefficient matrix with constant vector as last col. + * on return 'a' contains matrix inverse, last col contains + * solution vector. + */ +void mtrx_solv_gaussj (T **a, const size_t n); +void mtrx_solv_gaussj_v (T **a, T *v, const size_t n); +/** mtrx_solv_gaussj_inv returns mtrx_t containing inverse + solution. + * wrapper preserving original mtrx_t, allocating a copy before calling + * mtrx_solv_gaussj to create inverse of m and solution vector in newly + * allocted mtrx_t. user is responsible for freeing return. + */ +mtrx_t *mtrx_solv_gaussj_inv (const mtrx_t *m); +/** mtrx_get_sol_v return solution vector from inverse+solution mtrx_t. + * m must be a square matrix + constant vector as final column. returns + * allocated solution vector on success, NULL otherwise. + */ +T *mtrx_get_sol_v (const mtrx_t *m); + +/** qsort compare ascending for matrix */ +int mtrx_compare_rows_asc (const void *a, const void *b); + +#endif diff --git a/src/memrealloc.c b/src/memrealloc.c new file mode 100644 index 0000000..0b0459a --- /dev/null +++ b/src/memrealloc.c @@ -0,0 +1,51 @@ +#include "memrealloc.h" + +/** realloc 'ptr' of 'nelem' of 'psz' from 'oldsz' to 'newsz' of 'psz'. + * returns pointer to reallocated block of memory with new + * memory initialized to 0/NULL. return must be assigned to + * original pointer in caller. + */ +void *xrealloc_fixed (void *ptr, size_t psz, size_t oldsz, size_t newsz) +{ void *memptr = realloc ((char *)ptr, newsz * psz); + if (!memptr) { + perror ("realloc(): virtual memory exhausted."); + exit (EXIT_FAILURE); + } /* zero new memory (optional) */ + if (newsz > oldsz) + memset ((char *)memptr + oldsz * psz, 0, (newsz - oldsz) * psz); + return memptr; +} + +/** realloc 'ptr' of 'nelem' of 'psz' to 'nelem + inc' of 'psz'. + * returns pointer to reallocated block of memory with new + * memory initialized to 0/NULL. return must be assigned to + * original pointer in caller. + */ +void *xrealloc_inc (void *ptr, size_t psz, size_t *nelem, size_t inc) +{ void *memptr = realloc ((char *)ptr, (*nelem + inc) * psz); + if (!memptr) { + perror ("realloc(): virtual memory exhausted."); + exit (EXIT_FAILURE); + } /* zero new memory (optional) */ + memset ((char *)memptr + *nelem * psz, 0, inc * psz); + *nelem += inc; + return memptr; +} + +/** realloc 'ptr' of 'nelem' of 'psz' to 'nelem * 2' of 'psz'. + * returns pointer to reallocated block of memory with new + * memory initialized to 0/NULL. return must be assigned to + * original pointer in caller. + */ +void *xrealloc_x2 (void *ptr, size_t psz, size_t *nelem) +{ void *memptr = realloc ((char *)ptr, *nelem * 2 * psz); + if (!memptr) { + perror ("realloc(): virtual memory exhausted."); + exit (EXIT_FAILURE); + /* return NULL; */ + } /* zero new memory (optional) */ + memset ((char *)memptr + *nelem * psz, 0, *nelem * psz); + *nelem *= 2; + return memptr; +} + diff --git a/src/msolvtextview.c b/src/msolvtextview.c new file mode 100644 index 0000000..74059c6 --- /dev/null +++ b/src/msolvtextview.c @@ -0,0 +1,319 @@ +#include +#include "mtrx_t.h" + +#define PACKAGE "gtksolver" +#define VERSION "0.0.1" + +typedef struct app { + GtkWidget *text_view, /* text view */ + *btnsolv; /* solve button */ + GtkTextTag *tagblue, /* Deep sky blue tag */ + *tagroyal, /* Royal blue tag */ + *tagred; /* Red tag */ +} app_t; + +/* Solve... button callback */ +void btnsolv_activate (GtkWidget *widget, gpointer data) +{ + app_t *inst = data; + gint line; + size_t i; + gchar *buf, c; + gboolean havevalue = FALSE; + GtkTextIter start, end; + // GtkTextBuffer *buffer = GTK_TEXT_BUFFER(data); + GtkTextBuffer *buffer; + + buffer = gtk_text_view_get_buffer (GTK_TEXT_VIEW(inst->text_view)); + + /* get start and end iters for buffer */ + gtk_text_buffer_get_start_iter (buffer, &start); + gtk_text_buffer_get_end_iter (buffer, &end); + + /* advance start iter to beginning of first value */ + do { + c = gtk_text_iter_get_char (&start); + if (c == '.' || c == '-' || c == '+' || isdigit (c)) { + havevalue = TRUE; + break; + } + } while (gtk_text_iter_forward_char (&start) && + !gtk_text_iter_equal (&start, &end)); + + /* suppress incremental update of text_view */ + gtk_text_buffer_begin_user_action (buffer); + + /* get end iter for ouput */ + gtk_text_buffer_get_end_iter (buffer, &end); + + /* save current line to apply tag when done */ + line = gtk_text_iter_get_line (&end); + + if (havevalue) { + /* get text from buffer */ + buf = gtk_text_buffer_get_text (buffer, &start, &end, FALSE); + + /* fill matrix from values in buffer */ + mtrx_t *m = mtrx_read_alloc_buf (buf); + + /* solve system by gauss-jordan elimintation will full-pivoting */ + mtrx_solv_gaussj (m->mtrx, m->rows); /* m->mtrx now contains + * inverse + solution vector + */ + + /* output formatted solution vector */ + gtk_text_buffer_insert (buffer, &end, "\n\nSolution Vector:\n\n", -1); + for (i = 0; i < m->rows; i++) { /* output formatted solution vector */ + gchar *x; + x = g_strdup_printf (" x[%3zu] : % 11.7f\n", i, m->mtrx[i][m->rows]); + gtk_text_buffer_get_end_iter (buffer, &end); + gtk_text_buffer_insert (buffer, &end, x, -1); + g_free(x); + } + g_free (buf); + } + else + /* output format error - no system of equations found */ + gtk_text_buffer_insert (buffer, &end, + "\n\n => ERROR: No Numeric Value Found In Buffer\n", -1); + + /* get iters around solution vector title (lines +2 & 3 from save) */ + gtk_text_buffer_get_iter_at_line (buffer, &start, line+2); + gtk_text_buffer_get_iter_at_line (buffer, &end, line+3); + /* color output blue/red based on havevalue */ + gtk_text_buffer_apply_tag (buffer, + havevalue ? inst->tagblue: inst->tagred, + &start, &end); + + /* restore normal text_view behavior */ + gtk_text_buffer_end_user_action (buffer); + + + if(widget) {} +} + +/* Clear... button callback */ +void btnclear_activate (GtkWidget *widget, gpointer data) +{ + app_t *inst = data; + GtkTextBuffer *buffer; + buffer = gtk_text_view_get_buffer (GTK_TEXT_VIEW(inst->text_view)); + + gtk_text_buffer_set_text (buffer, "", -1); + + if (widget || data) {} +} + +/* Exit... button callback */ +void btnexit_activate (GtkWidget *widget, gpointer data) +{ + gtk_main_quit (); + + if (widget || data) {} +} + +static void settagtable (gpointer data) +{ + app_t *inst = data; + GtkTextBuffer *buffer; + buffer = gtk_text_view_get_buffer (GTK_TEXT_VIEW(inst->text_view)); + + /* tags to highlight button actions in intro + * X11 Color - Deep sky blue - "#00BFFF" + */ + inst->tagblue = gtk_text_buffer_create_tag (buffer, "blue_foreground", + "foreground", "#00BFFF", NULL); + /* X11 Color - Royal blue - "#4169E1" */ + inst->tagroyal = gtk_text_buffer_create_tag (buffer, "royal_foreground", + "foreground", "#4169E1", NULL); + /* X11 Color - Red */ + inst->tagred = gtk_text_buffer_create_tag (buffer, "red_foreground", + "foreground", "red", NULL); +} + +static void intro (gpointer data) +{ + app_t *inst = data; + GtkTextIter iter, end; + GtkTextTag *tagblue = inst->tagblue; + GtkTextBuffer *buffer; + buffer = gtk_text_view_get_buffer (GTK_TEXT_VIEW(inst->text_view)); + // GtkTextTag *tag/*, *tagg*/; + + gtk_text_buffer_get_iter_at_offset (buffer, &iter, 0); + // gtk_text_buffer_insert (buffer, &iter, "Linear System of Equations Solver\n\n", -1); + gtk_text_iter_forward_lines (&iter, 1); + gtk_text_buffer_insert (buffer, &iter, "\nClick [Clear] and enter Coefficient Matrix with Constant Vector as Last Column\n\n", -1); + gtk_text_iter_forward_lines (&iter, 1); + gtk_text_buffer_insert (buffer, &iter, "Example (3 x 4):\n\n" + " 3.0 2.0 -4.0 3.0\n" + " 2.0 3.0 3.0 15.0\n" + " 5.0 -3.0 1.0 14.0\n\nor\n\n" + "3,2,-4,3\n" + "2,3,3,15\n" + "5,-3,1,14\n\n", -1); + gtk_text_iter_forward_lines (&iter, 1); + gtk_text_buffer_insert (buffer, &iter, "Then click [Solve...]\n\n" + "Solution Vector:\n\n" + " x[ 0] : 3.0000000\n" + " x[ 1] : 1.0000000\n" + " x[ 2] : 2.0000000\n", -1); + + /* highlight [Clear...] */ + gtk_text_buffer_get_iter_at_line_offset (buffer, &iter, 1, 6); + gtk_text_buffer_get_iter_at_line_offset (buffer, &end, 1, 13); + gtk_text_buffer_apply_tag (buffer, tagblue, &iter, &end); + + /* highlight [Solve...] */ + gtk_text_buffer_get_iter_at_line_offset (buffer, &iter, 15, 11); + gtk_text_buffer_get_iter_at_line_offset (buffer, &end, 15, 21); + gtk_text_buffer_apply_tag (buffer, tagblue, &iter, &end); +} + +/* Help... button callback */ +void btnhelp_activate (GtkWidget *widget, gpointer data) +{ + app_t *inst = data; + GtkTextBuffer *buffer; + buffer = gtk_text_view_get_buffer (GTK_TEXT_VIEW(inst->text_view)); + + gtk_text_buffer_set_text (buffer, "", -1); + intro (data); + + if (widget || data) {} +} + +static void activate (app_t *inst) +{ + /* Declare variables */ + GtkWidget *window; + // GtkWidget *text_view; + GtkWidget *scrolled_window; + GtkWidget *vbox; + GtkWidget *hbox; + // GtkWidget *btnsolv, *btnclear, *btnexit; + GtkWidget *btnclear, *btnexit, *btnhelp; + GtkWidget *toplabel; + + GtkTextBuffer *buffer; + PangoFontDescription *font_desc; + // GdkColor color; + // GtkTextTag *tag/*, *tagg*/; + + /* Create a window with a title, and a default size (win, w, h) */ + window = gtk_window_new(GTK_WINDOW_TOPLEVEL); + gtk_window_set_title (GTK_WINDOW (window), "Linear System Solver"); + gtk_window_set_default_size (GTK_WINDOW (window), 650, 550); + + /* create vbox and add as main container in window */ + vbox = gtk_vbox_new (FALSE, 0); + gtk_container_add (GTK_CONTAINER (window), vbox); + + /* top label for title */ + toplabel = gtk_label_new (NULL); + gtk_label_set_markup (GTK_LABEL (toplabel), + "Matrix Solver for Linear System of Equations"); + gtk_misc_set_padding (GTK_MISC(toplabel), 0, 10); + gtk_box_pack_start (GTK_BOX (vbox), toplabel, FALSE, FALSE, 0); + gtk_widget_show (toplabel); + + /* Create the scrolled window. Usually NULL is passed for both parameters so + * that it creates the horizontal/vertical adjustments automatically. Setting + * the scrollbar policy to automatic allows the scrollbars to only show up + * when needed. + */ + scrolled_window = gtk_scrolled_window_new (NULL, NULL); + gtk_scrolled_window_set_policy (GTK_SCROLLED_WINDOW (scrolled_window), + GTK_POLICY_AUTOMATIC, + GTK_POLICY_AUTOMATIC); + gtk_container_set_border_width (GTK_CONTAINER (scrolled_window), 5); + gtk_box_pack_start(GTK_BOX(vbox), scrolled_window, TRUE, TRUE, 0); + + /* The text buffer represents for textview */ + buffer = gtk_text_buffer_new (NULL); + + /* Create text_view to display editable text buffer. + * Line wrapping is set to break lines in between words. + */ + inst->text_view = gtk_text_view_new_with_buffer (buffer); + // gtk_text_view_set_wrap_mode (GTK_TEXT_VIEW (inst->text_view), GTK_WRAP_WORD); + gtk_text_view_set_wrap_mode (GTK_TEXT_VIEW (inst->text_view), GTK_WRAP_NONE); + + /* Set default font throughout the widget */ + font_desc = pango_font_description_from_string ("DejaVu Sans Mono 8"); + gtk_widget_modify_font (inst->text_view, font_desc); + pango_font_description_free (font_desc); + + /* Change left margin throughout the widget */ + gtk_text_view_set_left_margin (GTK_TEXT_VIEW (inst->text_view), 10); + + /* Initialize tag table for text_view buffer */ + settagtable (inst); + + /* Add text_view to scrolled_window */ + gtk_container_add (GTK_CONTAINER (scrolled_window), inst->text_view); + + /* Create hbox to hold horizontal row of buttons at bottom of window */ + hbox = gtk_hbox_new (FALSE, 0); + gtk_box_set_spacing (GTK_BOX (hbox), 5); + gtk_container_set_border_width (GTK_CONTAINER (hbox), 5); + + /* Pack the hbox at the end of the vbox */ + gtk_box_pack_end (GTK_BOX(vbox), hbox, FALSE, FALSE, 0); + + /* define the buttons and pack into hbox */ + inst->btnsolv = gtk_button_new_with_mnemonic ("_Solve..."); + gtk_widget_set_size_request (inst->btnsolv, 80, 24); + // gtk_widget_set_sensitive (btnsolv, app->optregex); + gtk_box_pack_start (GTK_BOX (hbox), inst->btnsolv, FALSE, FALSE, 0); + + btnclear = gtk_button_new_with_mnemonic ("_Clear"); + gtk_widget_set_size_request (btnclear, 80, 24); + gtk_box_pack_start (GTK_BOX (hbox), btnclear, FALSE, FALSE, 0); + + btnhelp = gtk_button_new_with_mnemonic ("_Help"); + gtk_widget_set_size_request (btnhelp, 80, 24); + // gtk_box_pack_end (GTK_BOX (hbox), btnhelp, FALSE, FALSE, 0); + + btnexit = gtk_button_new_with_mnemonic ("_Exit"); + gtk_widget_set_size_request (btnexit, 80, 24); + gtk_box_pack_end (GTK_BOX (hbox), btnexit, FALSE, FALSE, 0); + /* push help after exit */ + gtk_box_pack_end (GTK_BOX (hbox), btnhelp, FALSE, FALSE, 0); + + /* connect window close callback */ + g_signal_connect (window, "destroy", + G_CALLBACK (gtk_main_quit), NULL); + + /* connect button callbacks */ + g_signal_connect (inst->btnsolv, "clicked", + G_CALLBACK (btnsolv_activate), inst); + + g_signal_connect (btnclear, "clicked", + G_CALLBACK (btnclear_activate), inst); + + g_signal_connect (btnhelp, "clicked", + G_CALLBACK (btnhelp_activate), inst); + + g_signal_connect (btnexit, "clicked", + G_CALLBACK (btnexit_activate), inst); + + /* Show intro text */ + intro (inst); + + gtk_widget_show_all (window); +} + + + +int main (int argc, char **argv) +{ + app_t inst = { .text_view = NULL }; + gtk_init(&argc, &argv); + + activate (&inst); + + gtk_main(); + + return 0; +} \ No newline at end of file diff --git a/src/mtrx_t.c b/src/mtrx_t.c new file mode 100644 index 0000000..964b50d --- /dev/null +++ b/src/mtrx_t.c @@ -0,0 +1,1171 @@ +#include "mtrx_t.h" +#include "memrealloc.h" + +/** debug - print all pointer values for mtrx_t->mtrx & mtrx_t->mtrx[i] */ +void prnptrs (mtrx_t *m) +{ + size_t i; + + printf ("m->mtrx: %p\n", (void*)m->mtrx); + for (i = 0; i < m->rowmax; i++) + printf (" m->mtrx[%2zu]: %p\n", i, (void*)m->mtrx[i]); + putchar ('\n'); +} + +/** allocate/initialize struct T and (m) pointers, set rowmax/colmax */ +mtrx_t *mtrx_create_ptrs (const size_t m) +{ + mtrx_t *ms = calloc (1, sizeof *ms); /* calloc struct */ + if (!ms) { + fprintf (stderr, "%s() error: calloc-m.\n", __func__); + return NULL; + } + + ms->mtrx = calloc (m, sizeof *ms->mtrx); /* calloc pointers */ + if (!ms->mtrx) { /* validate allocation */ + fprintf (stderr, "%s() error: malloc m->mtrx failed.\n", + __func__); + free (ms); + return NULL; + } + + ms->rowmax = m; + ms->colmax = COLSZ; + + return ms; +} + +/** allocate vector of n elements */ +T *vect_calloc (const size_t n) +{ + T *vect = calloc (n, sizeof *vect); /* calloc vector */ + + if (!vect) { + perror ("calloc-vector"); + fprintf (stderr, "%s() error: calloc-vect failed.\n", + __func__); + return NULL; + } + + return vect; +} + +/** allocate matrix of (m x n) */ +T **mtrx_calloc (const size_t m, const size_t n) +{ + register size_t i = 0; + T **mtrx = calloc (m, sizeof *mtrx); /* calloc pointers */ + if (!mtrx) { + perror ("calloc-pointers"); + fprintf (stderr, "%s() error: calloc mtrx failed.\n", + __func__); + return NULL; + } + + /* allocate m rows of n elements */ + for (i = 0; i < m; i++) { + if (!(mtrx[i] = calloc (n, sizeof *mtrx[i]))) { + perror ("calloc-mtrx[i]"); + fprintf (stderr, "%s() error: calloc-mtrx[%zu] failed.\n", + __func__, i); + while (i--) /* free all rows allocated */ + free (mtrx[i]); + free (mtrx); /* free pointers */ + return NULL; + } + } + + return mtrx; +} + +/** allocate/initialize struct T and (m x n) matrix, set rowmax/colmax */ +mtrx_t *mtrx_create_fixed (const size_t m, const size_t n) +{ + mtrx_t *ms = calloc (1, sizeof *ms); /* calloc struct */ + if (!ms) { + fprintf (stderr, "%s() error: calloc-m.\n", __func__); + return NULL; + } + + if (!(ms->mtrx = mtrx_calloc (m, n))) { /* allocate/validate matrix */ + free (ms); + return NULL; + } + + ms->rowmax = m; /* set rowmax/colmax to current allocation size */ + ms->colmax = n; + + return ms; /* return matrix struct */ +} + +/** realloc matrix pointers to twice current */ +T **mtrx_realloc_ptrs (mtrx_t *m) +{ + if (!m) { + fprintf (stderr, "%s() error: struct parameter NULL.\n", + __func__); + return NULL; + } + /* realloc m->mtrx to 2 x m->rowmax, zero new ptrs, set m->rowmax */ + m->mtrx = xrealloc_x2 (m->mtrx, sizeof *m->mtrx, &m->rowmax); + + return m->mtrx; +} + +/** realloc matrix_fixed, which = 0 realloc ptrs, 1 column storage. + * custom realloc used as no need to zero new pointers (they will be + * allocated) and row realloc simply saves the function call overhead + * for each row realloc. + */ +T **mtrx_realloc_fixed (mtrx_t *m, const int which) +{ + register size_t i; + + if (!m) { + fprintf (stderr, "%s() error: struct parameter NULL.\n", + __func__); + return NULL; + } + if (which == ROWOP) { /* reallocate pointers & new rows */ + void *tmp = realloc (m->mtrx, 2 * m->rowmax * sizeof *m->mtrx); + if (!tmp) { + fprintf (stderr, "%s() error: ptr realloc failed.\n", + __func__); + exit (EXIT_FAILURE); + } + m->mtrx = tmp; /* assign new block to m->mtrx */ + for (i = 0; i < m->rowmax; i++) { /* allocate rows */ + m->mtrx[i + m->rowmax] = calloc (m->colmax, sizeof **m->mtrx); + if (!m->mtrx[i + m->rowmax]) { + fprintf (stderr, "%s() error: ptr realloc failed.\n", + __func__); + exit (EXIT_FAILURE); + } + } + m->rowmax *= 2; /* incriment m->rowmax */ + return m->mtrx; + } + + /* reallocate row storage for all allocated rows */ + for (i = 0; i < m->rowmax; i++) { + void *tmp = realloc (m->mtrx[i], + 2 * m->colmax * sizeof *m->mtrx[i]); + if (!tmp) { + fprintf (stderr, "%s() error: realloc failed row[%zu].\n", + __func__, i); + exit (EXIT_FAILURE); + } + m->mtrx[i] = tmp; /* assign new block to (m->mtrx)[i] */ + memset (m->mtrx[i] + m->colmax, 0, /* zero new memory */ + m->colmax * sizeof *m->mtrx[i]); + } + m->colmax *= 2; /* increment colmax only if all succeed */ + + return m->mtrx; +} + +/** parse_dbl_array from any buffer regardless of format. + * uses strtod to parse all values from buf regardless of intervening + * characters. array initially sized to *nelem, reallocates as needed, + * a final realloc sizes the array to the number of elements, *nelem + * updated to final number of elements. returns allocated array on + * success, NULL otherwise. user responsible for freeing non-NULL return. + */ +T *parse_dbl_array (char *buf, size_t *nelem) +{ + T *array = NULL; + char *nptr = buf, *endptr; + size_t n = 0; + + if (!*nelem) + *nelem = COLSZ; + + if (!(array = calloc (*nelem, sizeof *array))) { + perror ("calloc-array"); + exit (EXIT_FAILURE); + } + + for (;;) { + errno = 0; /* reset errno before each conversion */ + if (n == *nelem) /* check if realloc required */ + array = xrealloc_x2 (array, sizeof *array, nelem); + /* skip any non-digit characters */ + while (*nptr && ((*nptr != '-' && *nptr != '+' && *nptr != '.' && + (*nptr < '0' || '9' < *nptr)) || + ((*nptr == '-' || *nptr == '+' || *nptr == '.') && + (*(nptr+1) < '0' || '9' < *(nptr+1))))) + nptr++; + + if (!*nptr) /* check if at end of buf */ + break; + array[n] = strtod (nptr, &endptr); /* call strtod */ + /* validate strod conversion */ + if (nptr == endptr) { /* no digits converted */ + fputs ("error: no digits converted.\n", stderr); + break; + } + else if (errno) { /* error in conversion under/overflow */ + fputs ("error: in conversion.\n", stderr); + break; + } + n++; /* increment element count */ + + /* skip delimiters/move pointer to next digit + * (avoids unnecessary realloc if no values remain) + */ + while (*endptr && *endptr != '-' && *endptr != '+' && *endptr != '.' + && (*endptr < '0' || *endptr > '9')) + endptr++; + if (!*endptr) /* break if end of string */ + break; + nptr = endptr; /* update pointer to end pointer */ + } + if (!n) { /* validate elements stored */ + free (array); + return NULL; + } + if (n < *nelem) { /* realloc to fit nrows, set rowmax */ + array = xrealloc_fixed (array, sizeof *array, *nelem, n); + *nelem = n; /* update rowmax */ + } + + return array; +} + +/** read (m x n) matrix from file stream into dynamically sized T *mtrx. + * pointers are reallocated x2 as needed, row storage is reallocated + * x2 as required to store all values in 1st row of data, then row storage + * is reallocate to shrink row allocation to exact number of columns and + * T->colmax is set. following read of all data, pointer storage is reduced + * to T->nrows and T->rowmax is set. + */ +mtrx_t *mtrx_read_alloc (FILE *fp) +{ + char buf[MAXC]; /* buffer for line */ + mtrx_t *m = mtrx_create_ptrs (ROWSZ); /* allocate ROWSIZE pointers */ + if (!m) /* validate */ + return NULL; + + while (fgets (buf, MAXC, fp)) { /* read each line */ + size_t col = 0; + + if (m->rows == m->rowmax) /* check if realloc needed */ + mtrx_realloc_ptrs (m); + + /* create array of doubles from buf */ + m->mtrx[m->rows] = parse_dbl_array (buf, &m->colmax); + if (!m->mtrx[m->rows]) { /* validate */ + fprintf (stderr, "%s() error: no values for mtrx[%zu][%zu].\n", + __func__, m->rows, m->cols); + exit (EXIT_FAILURE); + } + col = m->colmax; + + if (!m->rows) /* if first row */ + m->cols = col; /* set number of columns */ + else if (col != m->cols) /* check all rows have cols values */ + fprintf (stderr, "%s() error: column mismatch row[%zu]\n", + __func__, m->rows); + m->rows++; /* increment row count */ + } + + if (!m->rows) { + free (m->mtrx); + free (m); + return (m = NULL); + } + + if (m->rows < m->rowmax) { /* realloc to fit nrows, set rowmax */ + m->mtrx = xrealloc_fixed (m->mtrx, sizeof *m->mtrx, m->rowmax, m->rows); + m->rowmax = m->rows; /* update rowmax */ + } + + return m; /* return filled and exactly sized matrix struct */ +} + +mtrx_t *mtrx_read_alloc_buf (char *textviewbuf) +{ + char buf[MAXC]; /* buffer for line */ + char *p = textviewbuf; + size_t off = 0; + int used; + mtrx_t *m = mtrx_create_ptrs (ROWSZ); /* allocate ROWSIZE pointers */ + if (!m) /* validate */ + return NULL; + + while (sscanf (p + off, "%[^\n]%n", buf, &used) == 1) { /* parse each line */ + size_t col = 0; + + while (isspace (p[off + used])) + used++; + off += used; + + if (m->rows == m->rowmax) /* check if realloc needed */ + mtrx_realloc_ptrs (m); + + /* create array of doubles from buf */ + m->mtrx[m->rows] = parse_dbl_array (buf, &m->colmax); + if (!m->mtrx[m->rows]) { /* validate */ + fprintf (stderr, "%s() error: no values for mtrx[%zu][%zu].\n", + __func__, m->rows, m->cols); + exit (EXIT_FAILURE); + } + col = m->colmax; + + if (!m->rows) /* if first row */ + m->cols = col; /* set number of columns */ + else if (col != m->cols) /* check all rows have cols values */ + fprintf (stderr, "%s() error: column mismatch row[%zu]\n", + __func__, m->rows); + m->rows++; /* increment row count */ + } + + if (!m->rows) { + free (m->mtrx); + free (m); + return (m = NULL); + } + + if (m->rows < m->rowmax) { /* realloc to fit nrows, set rowmax */ + m->mtrx = xrealloc_fixed (m->mtrx, sizeof *m->mtrx, m->rowmax, m->rows); + m->rowmax = m->rows; /* update rowmax */ + } + + return m; /* return filled and exactly sized matrix struct */ +} + +/** matr_read_fixed, alloc/read (m x n) matrix from file stream. + * reallocate as required, but should not be required if correct initial + * size given. matrix is resized to exact number of rows and pointers + * before return. + */ +mtrx_t *mtrx_read_fixed (FILE *fp, const size_t rows, const size_t cols) +{ + char buf[MAXC]; + mtrx_t *m = mtrx_create_fixed (rows, cols); /* create m/alloc m->mtrx */ + if (!m) + return NULL; + + while (fgets (buf, MAXC, fp)) { /* read each line */ + char *nptr = buf, *endptr = buf; + size_t col = 0; + + if (m->rows == m->rowmax) /* check ptr realloc req'd */ + mtrx_realloc_fixed (m, ROWOP); /* custom realloc ptrs & rows */ + + /* NOTE cannot use parse_dbl_array () due to mtrx_realloc_fixed () + * (added skip of intervening characters) TODO - refactor. + */ + for (;;) { /* loop extracting all values in line */ + T v; /* value of type T */ + errno = 0; /* reset errno before conversion */ + /* skip any non-digit characters */ + while (*nptr && ((*nptr != '-' && *nptr != '+' && *nptr != '.' && + (*nptr < '0' || '9' < *nptr)) || + ((*nptr == '-' || *nptr == '+' || *nptr == '.') && + (*(nptr+1) < '0' || '9' < *(nptr+1))))) + nptr++; + + if (!*nptr) /* check if at end of buf */ + break; + v = strtod (nptr, &endptr); /* convert with strtod */ + if (nptr == endptr) { /* validate digits converted */ + fprintf (stderr, "%s() error: no digits for mtrx[%zu][%zu].\n", + __func__, m->rows, m->cols); + return m; + } + else if (errno) { /* check over/underflow in conversion */ + fprintf (stderr, "%s() error: failed conversion mtrx[%zu][%zu].\n", + __func__, m->rows, m->cols); + return m; + } + if (col == m->colmax) /* check row realloc req'd */ + mtrx_realloc_fixed (m, COLOP); /* custom row realloc */ + + m->mtrx[m->rows][col++] = v; /* assing value to row/col */ + + /* skip delimiters/move pointer to next digit + * (avoids unnecessary realloc if no values remain) + */ + while (*endptr && *endptr != '-' && *endptr != '+' && *endptr != '.' + && (*endptr < '0' || *endptr > '9')) + endptr++; + if (!*endptr) /* break if end of string */ + break; + nptr = endptr; /* update nptr from endptr */ + } + if (!m->rows) { /* set m->cols from 1st row */ + m->cols = col; + if (m->cols < m->colmax) { /* realloc to fit, set colmax */ + m->mtrx[m->rows] = xrealloc_fixed (m->mtrx[m->rows], + sizeof **m->mtrx, m->colmax, m->cols); + m->colmax = m->cols; /* update colmax */ + } + } + else if (col != m->cols) /* verify subsequent rows match */ + fprintf (stderr, "%s() error: column mismatch row[%zu]\n", + __func__, m->rows); + m->rows++; /* increment row count */ + } + if (m->rows < m->rowmax) { /* realloc to fit nrows, set rowmax */ + register size_t i; + for (i = m->rows; i < m->rowmax; i++) + free (m->mtrx[i]); + m->mtrx = xrealloc_fixed (m->mtrx, sizeof *m->mtrx, m->rowmax, m->rows); + m->rowmax = m->rows; /* update rowmax */ + } + + return m; +} + +/** print a (m x n) matrix */ +void mtrx_prn (const mtrx_t *m, int width) +{ + register size_t i, j; + + for (i = 0; i < m->rows; i++) { + for (j = 0; j < m->cols; j++) + printf (" % *g", width, m->mtrx[i][j]); + putchar ('\n'); + } +} + +/** print a (m x m+1) system of linear equations with separtor */ +void mtrx_sys_prn (const mtrx_t *m, int width) +{ + register size_t i, j; + + if (m->rows >= m->cols) { + fprintf (stderr, "%s() error - not a linear system matrix, " + "(rows >= cols).\n", __func__); + return; + } + + for (i = 0; i < m->rows; i++) { + for (j = 0; j < m->cols - 1; j++) + printf (" %*g", width, m->mtrx[i][j]); + printf (" | %*g\n", width, m->mtrx[i][m->cols - 1]); + } +} + +/** print a (m x m) matrix from (m x n) with n > m */ +void mtrx_prn_sq (const mtrx_t *m, int width) +{ + register size_t i, j; + + if (m->rows > m->cols) { + fprintf (stderr, "%s() error - square of m->rows not possible, " + "(rows > cols).\n", __func__); + return; + } + + for (i = 0; i < m->rows; i++) { + for (j = 0; j < m->rows; j++) + printf (" % *.3f", width, m->mtrx[i][j]); + putchar ('\n'); + } +} + +/** copy a struct matrix */ +mtrx_t *mtrx_copy (const mtrx_t *ma) +{ + register size_t i, j; + + mtrx_t *result = mtrx_create_fixed (ma->rows, ma->cols); + if (!result) /* allocate/validate stuct & (m x n) */ + return NULL; + + result->rowmax = result->rows = ma->rows; /* set retuls rows/cols */ + result->colmax = result->cols = ma->cols; + + /* copy mtrx contents */ + for (i = 0; i < ma->rows; i++) + for (j = 0; j < ma->cols; j++) + result->mtrx[i][j] = ma->mtrx[i][j]; + + return result; +} + +/** free a (m x n) matrix and struct for m->nrows pointers */ +void mtrx_free (mtrx_t *m) +{ + register size_t i; + + for (i = 0; i < m->rows; i++) + free (m->mtrx[i]); + free (m->mtrx); + free (m); +} + +void arr2d_free (T **m, const size_t n) +{ + register size_t i; + + for (i = 0; i < n; i++) + free (m[i]); + free (m); +} + +void arr_prn (T * const *a, const size_t m, const size_t n, const int wdth) +{ + register size_t i, j; + + for (i = 0; i < m; i++) { + for (j = 0; j < n; j++) + printf (" % *g", wdth, a[i][j]); + putchar ('\n'); + } +} + +void v_prn (const T *v, const size_t n, const int wdth) +{ + register size_t i; + + for (i = 0; i < n; i++) + printf (" % *g", wdth, v[i]); + + putchar ('\n'); +} + +/* check all elements of v are zero (within T_FP_MIN/MAX) */ +int v_is_zero_fp (const T *v, const size_t n) +{ + int is_zero = 1; + register size_t i; + + for (i = 0; i < n; i++) + if (v[i] < T_FP_MIN || T_FP_MAX < v[i]) { + is_zero = 0; + break; + } + + return is_zero; +} + +/** helper funcitons to incorporate in wrappers handling mtrx_t */ +/* copy two (m x n) matricies */ +T **m_copy (T * const *m_a, const size_t m, const size_t n) +{ + T **result = mtrx_calloc (m, n); + register size_t i, j; + + for (i = 0; i < m; i++) + for (j = 0; j < n; j++) + result [i][j] = m_a [i][j]; + + return result; +} + +/* add two (m x n) matricies together, return new result */ +T **m_add (T * const *m_a, T * const *m_b, const size_t m, const size_t n) +{ + T **result = mtrx_calloc (m, n); + register size_t i, j; + + for (i = 0; i < m; i++) + for (j = 0; j < n; j++) + result [i][j] = m_a [i][j] + m_b [i][j]; + + return result; +} + +/* subtract two (m x n) matricies together, return new result */ +T **m_sub (T * const *m_a, T * const *m_b, const size_t m, const size_t n) +{ + T **result = mtrx_calloc (m, n); + register size_t i, j; + + for (i = 0; i < m; i++) + for (j = 0; j < n; j++) + result [i][j] = m_a [i][j] - m_b [i][j]; + + return result; +} + +/* if _a is (m x n) and _b is (n x p), the product is an (m x p) matrix */ +T **m_mult (T * const *m_a, T * const *m_b, + const size_t m, const size_t n, const size_t p) +{ + T **result = mtrx_calloc (m, p); + register size_t i, j, k; + + for (i = 0; i < m; i++) + for (j = 0; j < p; j++) + for (k = 0; k < n; k++) + result [i][j] += m_a [i][k] * m_b [k][j]; + + return result; +} + +/* transpose matrix, dynamically allocate result */ +T **m_trans (T * const *m_a, size_t m, size_t n) +{ + T **result = mtrx_calloc (n, m); + register size_t i, j; + + for (i = 0; i < m; i++) + for (j = 0; j < n; j++) + result [j][i] = m_a [i][j]; + + return result; +} + +/* Transpose square matrix in place */ +void m_trans_sq (T **m_a, size_t n) +{ + register size_t i, j; + T tmp; + + for (i = 1; i < n; i++) { + for (j = 0; j < i; j++) { + tmp = m_a [i][j]; + m_a[i][j] = m_a[j][i]; + m_a[j][i] = tmp; + } + } +} + +/* cofactor matrix for (3 x 3) from matrix of minors */ +T **m_cofx (T * const *m, size_t n) +{ + T **cofx = mtrx_calloc (n, n); + register size_t i, j; + + /* fill cofactor from matrix of minors */ + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) { + int cof_1 = (i + 1) % n; + int cof_2 = (j + 1) % n; + int cof_3 = (i + 2) % n; + int cof_4 = (j + 2) % n; + + cofx[i][j]= (m[cof_1][cof_2] * m[cof_3][cof_4]) - + (m[cof_1][cof_4] * m[cof_3][cof_2]); + } + + return cofx; +} + +/* determinate for (3 x 3) from matrix and cofactor matrix */ +T m_det (T * const *m, T * const *cofx, const size_t n) +{ + T det = 0; + register size_t i; + + for (i = 0; i < n; i++) + det += m[i][0] * cofx[i][0]; + + return det; +} + +/* consistency check vector for (3 x 3) from adjugate and determinate */ +T *m_chk (T * const *adj, const T *v, const size_t n) +{ + T *chk = calloc (n, sizeof *chk); + register size_t i, j; + + if (!chk) { /* allocate/validate consistency check vector */ + fprintf (stderr, "%s() error: virtual memory exhausted.\n", + __func__); + return NULL; + } + + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) + chk[i] += adj[i][j] * v[j]; + + return chk; +} + +/* inverse matrix from determinate */ +T **m_inv (T * const *m, const T det, const size_t n) +{ + T **inv = mtrx_calloc (n, n); + register size_t i, j; + + if (!inv) + return NULL; + + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) + inv[i][j] = m[i][j] / det; + + return inv; +} + +/** add two struct matricies */ +mtrx_t *mtrx_add (const mtrx_t *ma, const mtrx_t *mb) +{ + mtrx_t *result = NULL; + + if (ma->rows != mb->rows || ma->cols != mb->cols) { + fprintf (stderr, "%s() error: unequal dimensions, (%zu x %zu) != " + "(%zu x %zu)\n", __func__, ma->rows, ma->cols, + mb->rows, mb->cols); + return NULL; + } + + /* alloccate/validate sturct */ + if (!(result = calloc (1, sizeof *result))) { + fprintf (stderr, "%s() error: calloc-result.\n", __func__); + return NULL; + } + + result->rowmax = result->rows = ma->rows; /* set retuls rows/cols */ + result->colmax = result->cols = ma->cols; + + /* add ma->mtrx and mb->mtrx assigning newly allocated result */ + result->mtrx = m_add (ma->mtrx, mb->mtrx, ma->rows, ma->cols); + + return result; +} + +/** subtract two struct matricies */ +mtrx_t *mtrx_sub (const mtrx_t *ma, const mtrx_t *mb) +{ + mtrx_t *result = NULL; + + if (ma->rows != mb->rows || ma->cols != mb->cols) { + fprintf (stderr, "%s() error: unequal dimensions, (%zu x %zu) != " + "(%zu x %zu)\n", __func__, ma->rows, ma->cols, + mb->rows, mb->cols); + return NULL; + } + + /* alloccate/validate sturct */ + if (!(result = calloc (1, sizeof *result))) { + fprintf (stderr, "%s() error: calloc-result.\n", __func__); + return NULL; + } + + result->rowmax = result->rows = ma->rows; /* set retuls rows/cols */ + result->colmax = result->cols = ma->cols; + + /* add ma->mtrx and mb->mtrx assigning newly allocated result */ + result->mtrx = m_sub (ma->mtrx, mb->mtrx, ma->rows, ma->cols); + + return result; +} + +/** multiply two struct matricies */ +mtrx_t *mtrx_mult (const mtrx_t *ma, const mtrx_t *mb) +{ + mtrx_t *result = NULL; + + if (ma->cols != mb->rows) { + fprintf (stderr, "%s() error: unequal dimensions, (%zu x %zu) != " + "(%zu x %zu)\n [not (m x n) * (n x p) => (m x p)]\n", + __func__, ma->rows, ma->cols, mb->rows, mb->cols); + return NULL; + } + + /* alloccate/validate sturct */ + if (!(result = calloc (1, sizeof *result))) { + fprintf (stderr, "%s() error: calloc-result.\n", __func__); + return NULL; + } + + result->rowmax = result->rows = mb->rows; /* set retuls rows/cols */ + result->colmax = result->cols = mb->cols; + + /* add ma->mtrx and mb->mtrx assigning newly allocated result */ + result->mtrx = m_mult (ma->mtrx, mb->mtrx, ma->rows, ma->cols, mb->cols); + + return result; +} + +/** transpose matrix */ +mtrx_t *mtrx_trans (const mtrx_t *m) +{ + mtrx_t *result = NULL; + + /* alloccate/validate sturct */ + if (!(result = calloc (1, sizeof *result))) { + fprintf (stderr, "%s() error: calloc-result.\n", __func__); + return NULL; + } + + result->rowmax = result->rows = m->cols; /* set retuls rows/cols */ + result->colmax = result->cols = m->rows; + + /* add ma->mtrx and mb->mtrx assigning newly allocated result */ + result->mtrx = m_trans (m->mtrx, m->rows, m->cols); + + return result; +} + +/** solve system of equations - FIXME cofactor for N not 3x3 + * (just use src-c/tmp/gauss_jordan_chk_10x10_mod.c) and change + * to accept mtrx_t *m + solution vect. + */ +T *mtrx_solv (const mtrx_t *m, const T *v) +{ + T **cofx = NULL, /* cofactor matrix (adjugate as transposed) */ + **inv = NULL; /* inverse of adjugate matrix */ + T *chk = NULL, /* consistency check vector */ + *sol = NULL; /* solution vector */ + T det = 0; /* determinate */ + register size_t i, j; /* loop vars */ + + if (m->rows != m->cols) { /* validate m is a square matrix */ + fprintf (stderr, "%s() error: maxtix not square (%zu x %zu).\n", + __func__, m->rows, m->cols); + return NULL; + } + + /* allocate/validate solution vector */ + if (!(sol = calloc (m->rows, sizeof *sol))) { + fprintf (stderr, "%s() error: memory exhausted 'sol'.\n", __func__); + return NULL; + } + + /* form cofactor matrix */ + if (!(cofx = m_cofx (m->mtrx, m->rows))) + return NULL; + + /* find determinate */ + det = m_det (m->mtrx, cofx, m->rows); +#ifdef DEBUG + puts ("\ncoefficient matrix:\n"); + arr_prn (m->mtrx, m->rows, m->rows, 5); + + puts ("\ncofactor matrix:\n"); + arr_prn (cofx, m->rows, m->rows, 5); + + printf ("\ndeterminant: %.2f\n", det); +#endif + /* form adjugate with in-place transpose of cofactor */ + m_trans_sq (cofx, m->rows); +#ifdef DEBUG + puts ("\nadjugate matrix:\n"); + arr_prn (cofx, m->rows, m->rows, 5); +#endif + /* form consistency check vector */ + if (!(chk = m_chk (cofx, v, m->rows))) { + arr2d_free (cofx, m->rows); + return NULL; + } +#ifdef DEBUG + puts ("\nconsistency check vector:\n"); + v_prn (chk, m->rows, 5); + putchar ('\n'); +#endif + if (det == 0) { /* eliminate infinite and no solution cases */ + if (v_is_zero_fp (chk, m->rows)) /* consistency vect zero */ + printf ("The system is consistent - infinite solutions.\n"); + else + printf ("The system has no solution.\n"); + + arr2d_free (cofx, m->rows); + free (chk); + + return NULL; + } + + /* form inverse of adjugate matrix */ + if (!(inv = m_inv (cofx, det, m->rows))) { + arr2d_free (cofx, m->rows); + free (chk); + return NULL; + } + +#ifdef DEBUG + puts ("\ninverse of adjugate;\n"); + arr_prn (inv, m->rows, m->rows, 5); +#endif + /* compute solution vector by inverse */ + for (i = 0; i < m->rows; i++) + for (j = 0; j < m->rows; j++) + sol[i] += inv[i][j] * v[j]; + + arr2d_free (inv, m->rows); + arr2d_free (cofx, m->rows); + free (chk); + + /* check if solution vector is zero - trivial solution */ + if (v_is_zero_fp (sol, m->rows)) { + printf ("The system has a trivial solution.\n"); + free (sol); + return NULL; + } + + return sol; /* return unique solution vector */ +} + +/** mtrx_solv_cmb - solves system of eq where m contains solution vect. + * solves system where m contains the solution vector as the last column + * in the form {m} = {m'}[v] where {m'} is the coefficient matrix and [v] + * the solution vector. it is a wrapper for mtrx_solv() above. returns + * a vector containing the unique solution or NULL if {m'} is singular + * or the solutions are infinite or trivial. + */ +T *mtrx_solv_cmb (mtrx_t *m) +{ + mtrx_t *mc = NULL; /* matrix of coefficients */ + T *v = NULL, /* vector of constants */ + *sol = NULL; /* unique solution vector */ + register size_t i; + + /* create matrix struct for mc */ + if (!(mc = calloc (1, sizeof *mc))) + return NULL; + + /* create/allocaate/copy coefficient maxtrix */ + if (!(mc->mtrx = m_copy (m->mtrx, m->rows, m->cols - 1))) + return NULL; + + mc->rows = m->rows; /* set rows/cols */ + mc->cols = m->cols - 1; + + if (!(v = calloc (m->rows, sizeof *v))) { + fprintf (stderr, "%s() error: calloc-v\n", __func__); + return NULL; + } + + for (i = 0; i < m->rows; i++) + v[i] = m->mtrx[i][m->cols-1]; + + sol = mtrx_solv (mc, v); + + mtrx_free (mc); + free (v); + + return sol; +} + +static void SWAP (T *a, T *b) +{ + T tmp = *a; + *a = *b; + *b = tmp; +} + +/** Guass-Jordan elimination with full pivoting. + * 'a' is coefficient matrix with constant vector as last col. + * on return 'a' contains matrix inverse, last col contains + * solution vector. + */ +void mtrx_solv_gaussj (T **a, const size_t n) +{ /* bookkeeping arrays for pivot */ + int *indxc = calloc (n, sizeof *indxc), + *indxr = calloc (n, sizeof *indxr), + *ipiv = calloc (n, sizeof *ipiv); + T big, dum, pivinv; + size_t i, icol = 0, irow = 0, j, k, l, ll; + + for (j = 0; j < n; j++) ipiv[j] = 0; + for (i = 0; i < n; i++) { + big = 0.0; + for (j = 0; j < n; j++) + if (ipiv[j] != 1) + for (k = 0; k < n; k++) { + if (ipiv[k] == 0) { + if (fabs (a[j][k]) >= big) { + big = fabs (a[j][k]); + irow = j; + icol = k; + } + } + } + ipiv[icol]++; + + if (irow != icol) /* transpose row/col */ + for (l = 0; l < n + 1; l++) SWAP(&a[irow][l], &a[icol][l]); + + indxr[i] = irow; + indxc[i] = icol; + + if (a[icol][icol] == 0.0) { + /* nrerror ("gaussj: Singular Matrix"); */ + fprintf (stderr, "guassj() error: singular matrix.\n"); + goto gaussjdone; + } + + pivinv = 1.0 / a[icol][icol]; + a[icol][icol] = 1.0; + + for (l = 0; l < n + 1; l++) a[icol][l] *= pivinv; + + for (ll = 0; ll < n; ll++) + if (ll != icol) { + dum = a[ll][icol]; + a[ll][icol] = 0.0; + for (l = 0; l < n + 1; l++) a[ll][l] -= a[icol][l] * dum; + } + } + + l = n; + while (l--) { + if (indxr[l] != indxc[l]) + for (k = 0; k < n; k++) + SWAP(&a[k][indxr[l]],&a[k][indxc[l]]); + } + gaussjdone:; + + free (ipiv); + free (indxr); + free (indxc); +} + +/** same taking a (n x n) and v. TODO preserve T * const * a by making a copy + * allocating for the inverse and returing the inverse while letting the + * solution vector be available to the caller through the updated v. + */ +void mtrx_solv_gaussj_v (T **a, T *v, const size_t n) +{ /* bookkeeping arrays for pivot */ + int *indxc = calloc (n, sizeof *indxc), + *indxr = calloc (n, sizeof *indxr), + *ipiv = calloc (n, sizeof *ipiv); + T big, dum, pivinv; + size_t i, icol = 0, irow = 0, j, k, l, ll; + + for (j = 0; j < n; j++) ipiv[j] = 0; + for (i = 0; i < n; i++) { + big = 0.0; + for (j = 0; j < n; j++) + if (ipiv[j] != 1) + for (k = 0; k < n; k++) { + if (ipiv[k] == 0) { + if (fabs (a[j][k]) >= big) { + big = fabs (a[j][k]); + irow = j; + icol = k; + } + } + } + ipiv[icol]++; + + if (irow != icol) { /* transpose row/col */ + for (l = 0; l < n; l++) SWAP(&a[irow][l], &a[icol][l]); + SWAP (&v[irow], &v[icol]); + } + indxr[i] = irow; + indxc[i] = icol; + + if (a[icol][icol] == 0.0) { + /* nrerror ("gaussj: Singular Matrix"); */ + fprintf (stderr, "guassj() error: singular matrix.\n"); + goto gaussjdone; + } + + pivinv = 1.0 / a[icol][icol]; + a[icol][icol] = 1.0; + + for (l = 0; l < n; l++) a[icol][l] *= pivinv; + v[icol] *= pivinv; + + for (ll = 0; ll < n; ll++) + if (ll != icol) { + dum = a[ll][icol]; + a[ll][icol] = 0.0; + for (l = 0; l < n; l++) a[ll][l] -= a[icol][l] * dum; + v[ll] -= v[icol] * dum; + } + } + + l = n; + while (l--) { + if (indxr[l] != indxc[l]) + for (k = 0; k < n; k++) + SWAP(&a[k][indxr[l]],&a[k][indxc[l]]); + } + gaussjdone:; + + free (ipiv); + free (indxr); + free (indxc); +} + +/** mtrx_solv_gaussj_inv returns mtrx_t containing inverse + solution. + * wrapper preserving original mtrx_t, allocating a copy before calling + * mtrx_solv_gaussj to create inverse of m and solution vector in newly + * allocted mtrx_t. user is responsible for freeing return. + */ +mtrx_t *mtrx_solv_gaussj_inv (const mtrx_t *m) +{ + mtrx_t *invsol = NULL; + + if (m->cols <= m->rows) { + fprintf (stderr, "%s() error: invalid mtrx_t size (cols <= rows)\n", + __func__); + return NULL; + } + + if (!(invsol = mtrx_copy (m))) + return NULL; + + mtrx_solv_gaussj (invsol->mtrx, invsol->rows); + + return invsol; +} + +/** mtrx_get_sol_v return solution vector from inverse+solution mtrx_t. + * m must be a square matrix + constant vector as final column. returns + * allocated solution vector on success, NULL otherwise. + */ +T *mtrx_get_sol_v (const mtrx_t *m) +{ + T *sol = NULL; + register size_t i; + + if (m->cols <= m->rows) { + fprintf (stderr, "%s() error: invalid mtrx_t size (cols <= rows)\n", + __func__); + return NULL; + } + + if (!(sol = calloc (m->rows, sizeof *sol))) { + fprintf (stderr, "%s() error: memory exhausted calloc-sol.\n", + __func__); + return NULL; + } + + for (i = 0; i < m->rows; i++) + sol[i] = m->mtrx[i][m->rows]; + + return sol; +} + +/** qsort compare ascending for matrix (1st value only) + * must write a wrapper to sort by rows by col values, + * and then sort by 1st value. + */ +int mtrx_compare_rows_asc (const void *a, const void *b) +{ + /* (a > b) - (a < b) */ + return (**(T * const *)a > **(T * const *)b) - + (**(T * const *)a < **(T * const *)b); +} + + +// void *tmp = realloc (m->mtrx[m->rows], +// 2 * m->colmax * sizeof **m->mtrx); +// if (!tmp) { /* validate */ +// perror ("realloc-m->mtrx[m->rows]"); +// exit (EXIT_FAILURE); +// } +// m->mtrx[m->rows] = tmp; /* assigne new block of memory */ +// memset (m->mtrx[m->rows] + m->colmax, 0, /* zero new memory */ +// m->colmax * sizeof *m->mtrx[m->rows]); +// m->colmax *= 2; /* increment colmax */ + +// void *tmp = realloc (m->mtrx[m->rows], +// m->cols * sizeof **m->mtrx); +// m->colmax = m->cols; /* update colmax */ +// if (!tmp) { /* validate allocation */ +// perror ("realloc-m->mtrx[m->rows]_to_fit"); +// exit (EXIT_FAILURE); +// } +// m->mtrx[m->rows] = tmp; /* assign new block of memory */ + +// mtrx_realloc_ptrs (T *m) +// void *tmp = realloc (m->mtrx, 2 * m->rowmax * sizeof *m->mtrx); +// if (!tmp) { +// fprintf (stderr, "%s() error: ptr realloc failed.\n", +// __func__); +// exit (EXIT_FAILURE); +// } +// m->mtrx = tmp; /* assign new block to m->mtrx */ +// memset (m->mtrx + m->rowmax, 0, /* zero new pointers */ +// m->rowmax * sizeof *m->mtrx); +// m->rowmax *= 2; /* incriment m->rowmax */