Changeset 9e269a0 in git for Singular/extra.cc
 Timestamp:
 Dec 22, 2010, 1:54:29 PM (13 years ago)
 Branches:
 (u'spielwiese', 'd0474371d8c5d8068ab70bfb42719c97936b18a6')
 Children:
 e2fd5d453ccf5a7982b273cf039c5c7f41f41080
 Parents:
 2fd30cb3d4e41a93238b92df65a30fb88c8276de
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/extra.cc
r2fd30cb r9e269a0 667 667 } 668 668 } 669 /*==================== lduDecomp ======================*/ 670 if(strcmp(sys_cmd, "lduDecomp")==0) 671 { 672 if ((h != NULL) && (h>Typ() == MATRIX_CMD) && (h>next == NULL)) 673 { 674 matrix aMat = (matrix)h>Data(); 675 matrix pMat; matrix lMat; matrix dMat; matrix uMat; 676 poly l; poly u; poly prodLU; 677 lduDecomp(aMat, pMat, lMat, dMat, uMat, l, u, prodLU); 678 lists L = (lists)omAllocBin(slists_bin); 679 L>Init(7); 680 L>m[0].rtyp = MATRIX_CMD; L>m[0].data=(void*)pMat; 681 L>m[1].rtyp = MATRIX_CMD; L>m[1].data=(void*)lMat; 682 L>m[2].rtyp = MATRIX_CMD; L>m[2].data=(void*)dMat; 683 L>m[3].rtyp = MATRIX_CMD; L>m[3].data=(void*)uMat; 684 L>m[4].rtyp = POLY_CMD; L>m[4].data=(void*)l; 685 L>m[5].rtyp = POLY_CMD; L>m[5].data=(void*)u; 686 L>m[6].rtyp = POLY_CMD; L>m[6].data=(void*)prodLU; 687 res>rtyp = LIST_CMD; 688 res>data = (char *)L; 689 return FALSE; 690 } 691 else 692 { 693 Werror( "expected argument list (int, int, poly, poly, poly, int)"); 694 return TRUE; 695 } 696 } 697 /*==================== lduSolve ======================*/ 698 if(strcmp(sys_cmd, "lduSolve")==0) 699 { 700 /* for solving a linear equation system A * x = b, via the 701 given LDUdecomposition of the matrix A; 702 There is one valid parametrisation: 703 1) exactly eight arguments P, L, D, U, l, u, lTimesU, b; 704 P, L, D, and U realise the LDUdecomposition of A, that is, 705 P * A = L * D^(1) * U, and P, L, D, and U satisfy the 706 properties decribed in method 'luSolveViaLDUDecomp' in 707 linearAlgebra.h; see there; 708 l, u, and lTimesU are as described in the same location; 709 b is the righthand side vector of the linear equation system; 710 The method will return a list of either 1 entry or three entries: 711 1) [0] if there is no solution to the system; 712 2) [1, x, H] if there is at least one solution; 713 x is any solution of the given linear system, 714 H is the matrix with column vectors spanning the homogeneous 715 solution space. 716 The method produces an error if matrix and vector sizes do not 717 fit. */ 718 if ((h == NULL)  (h>Typ() != MATRIX_CMD)  719 (h>next == NULL)  (h>next>Typ() != MATRIX_CMD)  720 (h>next>next == NULL)  (h>next>next>Typ() != MATRIX_CMD)  721 (h>next>next>next == NULL)  722 (h>next>next>next>Typ() != MATRIX_CMD)  723 (h>next>next>next>next == NULL)  724 (h>next>next>next>next>Typ() != POLY_CMD)  725 (h>next>next>next>next>next == NULL)  726 (h>next>next>next>next>next>Typ() != POLY_CMD)  727 (h>next>next>next>next>next>next == NULL)  728 (h>next>next>next>next>next>next>Typ() != POLY_CMD)  729 (h>next>next>next>next>next>next>next == NULL)  730 (h>next>next>next>next>next>next>next>Typ() 731 != MATRIX_CMD)  732 (h>next>next>next>next>next>next>next>next != NULL)) 733 { 734 Werror("expected input (matrix, matrix, matrix, matrix, %s", 735 "poly, poly, poly, matrix)"); 736 return TRUE; 737 } 738 matrix pMat = (matrix)h>Data(); 739 matrix lMat = (matrix)h>next>Data(); 740 matrix dMat = (matrix)h>next>next>Data(); 741 matrix uMat = (matrix)h>next>next>next>Data(); 742 poly l = (poly) h>next>next>next>next>Data(); 743 poly u = (poly) h>next>next>next>next>next>Data(); 744 poly lTimesU = (poly) h>next>next>next>next>next>next 745 >Data(); 746 matrix bVec = (matrix)h>next>next>next>next>next>next 747 >next>Data(); 748 matrix xVec; int solvable; matrix homogSolSpace; 749 if (pMat>rows() != pMat>cols()) 750 { 751 Werror("first matrix (%d x %d) is not quadratic", 752 pMat>rows(), pMat>cols()); 753 return TRUE; 754 } 755 if (lMat>rows() != lMat>cols()) 756 { 757 Werror("second matrix (%d x %d) is not quadratic", 758 lMat>rows(), lMat>cols()); 759 return TRUE; 760 } 761 if (dMat>rows() != dMat>cols()) 762 { 763 Werror("third matrix (%d x %d) is not quadratic", 764 dMat>rows(), dMat>cols()); 765 return TRUE; 766 } 767 if (dMat>cols() != uMat>rows()) 768 { 769 Werror("third matrix (%d x %d) and fourth matrix (%d x %d) %s", 770 "do not t", 771 dMat>rows(), dMat>cols(), uMat>rows(), uMat>cols()); 772 return TRUE; 773 } 774 if (uMat>rows() != bVec>rows()) 775 { 776 Werror("fourth matrix (%d x %d) and vector (%d x 1) do not fit", 777 uMat>rows(), uMat>cols(), bVec>rows()); 778 return TRUE; 779 } 780 solvable = luSolveViaLDUDecomp(pMat, lMat, dMat, uMat, l, u, lTimesU, 781 bVec, xVec, homogSolSpace); 782 783 /* build the return structure; a list with either one or 784 three entries */ 785 lists ll = (lists)omAllocBin(slists_bin); 786 if (solvable) 787 { 788 ll>Init(3); 789 ll>m[0].rtyp=INT_CMD; ll>m[0].data=(void *)solvable; 790 ll>m[1].rtyp=MATRIX_CMD; ll>m[1].data=(void *)xVec; 791 ll>m[2].rtyp=MATRIX_CMD; ll>m[2].data=(void *)homogSolSpace; 792 } 793 else 794 { 795 ll>Init(1); 796 ll>m[0].rtyp=INT_CMD; ll>m[0].data=(void *)solvable; 797 } 798 res>rtyp = LIST_CMD; 799 res>data=(char*)ll; 800 return FALSE; 801 } 669 802 /*==================== forking experiments ======================*/ 670 803 if(strcmp(sys_cmd, "waitforssilinks")==0)
Note: See TracChangeset
for help on using the changeset viewer.