Mercurial > octave

comparison liboctave/CSparse.cc @ 5203:dbeafbc0ff64

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-03-15 00:58:55 by jwe]
author jwe
date Tue, 15 Mar 2005 00:58:56 +0000
parents 57077d0ddc8e
children 90a9058de7e8
comparison
equal deleted inserted replaced
5202:49de4e624020 5203:dbeafbc0ff64
38 #include "boolSparse.h" 38 #include "boolSparse.h"
39 #include "dSparse.h" 39 #include "dSparse.h"
40 #include "oct-spparms.h" 40 #include "oct-spparms.h"
41 #include "SparseCmplxLU.h" 41 #include "SparseCmplxLU.h"
42 42
43 #ifdef HAVE_UMFPACK
43 // External UMFPACK functions in C 44 // External UMFPACK functions in C
44 extern "C" { 45 extern "C" {
45 #include "umfpack.h" 46 #include <umfpack/umfpack.h>
46 } 47 }
48 #endif
47 49
48 // Fortran functions we call. 50 // Fortran functions we call.
49 extern "C" 51 extern "C"
50 { 52 {
51 F77_RET_T 53 F77_RET_T
625 627
626 ComplexDET 628 ComplexDET
627 SparseComplexMatrix::determinant (int& err, double& rcond, int calc_cond) const 629 SparseComplexMatrix::determinant (int& err, double& rcond, int calc_cond) const
628 { 630 {
629 ComplexDET retval; 631 ComplexDET retval;
632 #ifdef HAVE_UMFPACK
630 633
631 int nr = rows (); 634 int nr = rows ();
632 int nc = cols (); 635 int nc = cols ();
633 636
634 if (nr == 0 || nc == 0 || nr != nc) 637 if (nr == 0 || nc == 0 || nr != nc)
740 else 743 else
741 retval = ComplexDET (d); 744 retval = ComplexDET (d);
742 } 745 }
743 } 746 }
744 } 747 }
748 #else
749 (*current_liboctave_error_handler) ("UMFPACK not installed");
750 #endif
745 751
746 return retval; 752 return retval;
747 } 753 }
748 754
749 ComplexMatrix 755 ComplexMatrix
4444 { 4450 {
4445 // The return values 4451 // The return values
4446 void *Numeric; 4452 void *Numeric;
4447 err = 0; 4453 err = 0;
4448 4454
4455 #ifdef HAVE_UMFPACK
4449 // Setup the control parameters 4456 // Setup the control parameters
4450 Control = Matrix (UMFPACK_CONTROL, 1); 4457 Control = Matrix (UMFPACK_CONTROL, 1);
4451 double *control = Control.fortran_vec (); 4458 double *control = Control.fortran_vec ();
4452 umfpack_zi_defaults (control); 4459 umfpack_zi_defaults (control);
4453 4460
4540 } 4547 }
4541 } 4548 }
4542 4549
4543 if (err != 0) 4550 if (err != 0)
4544 umfpack_zi_free_numeric (&Numeric); 4551 umfpack_zi_free_numeric (&Numeric);
4552 #else
4553 (*current_liboctave_error_handler) ("UMFPACK not installed");
4554 #endif
4545 4555
4546 return Numeric; 4556 return Numeric;
4547 } 4557 }
4548 4558
4549 ComplexMatrix 4559 ComplexMatrix
4578 typ = SparseType::Full; 4588 typ = SparseType::Full;
4579 } 4589 }
4580 4590
4581 if (typ == SparseType::Full) 4591 if (typ == SparseType::Full)
4582 { 4592 {
4593 #ifdef HAVE_UMFPACK
4583 Matrix Control, Info; 4594 Matrix Control, Info;
4584 void *Numeric = factorize (err, rcond, Control, Info, 4595 void *Numeric = factorize (err, rcond, Control, Info,
4585 sing_handler); 4596 sing_handler);
4586 4597
4587 if (err == 0) 4598 if (err == 0)
4592 double *control = Control.fortran_vec (); 4603 double *control = Control.fortran_vec ();
4593 double *info = Info.fortran_vec (); 4604 double *info = Info.fortran_vec ();
4594 const int *Ap = cidx (); 4605 const int *Ap = cidx ();
4595 const int *Ai = ridx (); 4606 const int *Ai = ridx ();
4596 const Complex *Ax = data (); 4607 const Complex *Ax = data ();
4608 #ifdef UMFPACK_SEPARATE_SPLIT
4597 const double *Bx = b.fortran_vec (); 4609 const double *Bx = b.fortran_vec ();
4598 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); 4610 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
4599 for (int i = 0; i < b_nr; i++) 4611 for (int i = 0; i < b_nr; i++)
4600 Bz[i] = 0.; 4612 Bz[i] = 0.;
4601 4613 #else
4614 OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr);
4615 #endif
4602 retval.resize (b_nr, b_nc); 4616 retval.resize (b_nr, b_nc);
4603 Complex *Xx = retval.fortran_vec (); 4617 Complex *Xx = retval.fortran_vec ();
4604 4618
4605 for (int j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) 4619 for (int j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
4606 { 4620 {
4621 #ifdef UMFPACK_SEPARATE_SPLIT
4607 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, 4622 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai,
4608 X_CAST (const double *, Ax), 4623 X_CAST (const double *, Ax),
4609 NULL, 4624 NULL,
4610 X_CAST (double *, &Xx[iidx]), 4625 X_CAST (double *, &Xx[iidx]),
4611 NULL, 4626 NULL,
4612 &Bx[iidx], Bz, Numeric, 4627 &Bx[iidx], Bz, Numeric,
4613 control, info); 4628 control, info);
4629 #else
4630 for (int i = 0; i < b_nr; i++)
4631 Bz[i] = b.elem (i, j);
4632
4633 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai,
4634 X_CAST (const double *, Ax),
4635 NULL,
4636 X_CAST (double *, &Xx[iidx]),
4637 NULL,
4638 X_CAST (const double *, Bz),
4639 NULL, Numeric,
4640 control, info);
4641 #endif
4642
4614 if (status < 0) 4643 if (status < 0)
4615 { 4644 {
4616 (*current_liboctave_error_handler) 4645 (*current_liboctave_error_handler)
4617 ("SparseComplexMatrix::solve solve failed"); 4646 ("SparseComplexMatrix::solve solve failed");
4618 4647
4645 4674
4646 umfpack_zi_report_info (control, info); 4675 umfpack_zi_report_info (control, info);
4647 4676
4648 umfpack_zi_free_numeric (&Numeric); 4677 umfpack_zi_free_numeric (&Numeric);
4649 } 4678 }
4679 #else
4680 (*current_liboctave_error_handler) ("UMFPACK not installed");
4681 #endif
4650 } 4682 }
4651 else if (typ != SparseType::Hermitian) 4683 else if (typ != SparseType::Hermitian)
4652 (*current_liboctave_error_handler) ("incorrect matrix type"); 4684 (*current_liboctave_error_handler) ("incorrect matrix type");
4653 } 4685 }
4654 4686
4687 typ = SparseType::Full; 4719 typ = SparseType::Full;
4688 } 4720 }
4689 4721
4690 if (typ == SparseType::Full) 4722 if (typ == SparseType::Full)
4691 { 4723 {
4724 #ifdef HAVE_UMFPACK
4692 Matrix Control, Info; 4725 Matrix Control, Info;
4693 void *Numeric = factorize (err, rcond, Control, Info, sing_handler); 4726 void *Numeric = factorize (err, rcond, Control, Info, sing_handler);
4694 4727
4695 if (err == 0) 4728 if (err == 0)
4696 { 4729 {
4701 double *info = Info.fortran_vec (); 4734 double *info = Info.fortran_vec ();
4702 const int *Ap = cidx (); 4735 const int *Ap = cidx ();
4703 const int *Ai = ridx (); 4736 const int *Ai = ridx ();
4704 const Complex *Ax = data (); 4737 const Complex *Ax = data ();
4705 4738
4739 #ifdef UMFPACK_SEPARATE_SPLIT
4706 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); 4740 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
4707 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); 4741 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
4708 for (int i = 0; i < b_nr; i++) 4742 for (int i = 0; i < b_nr; i++)
4709 Bz[i] = 0.; 4743 Bz[i] = 0.;
4744 #else
4745 OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr);
4746 #endif
4710 4747
4711 // Take a first guess that the number of non-zero terms 4748 // Take a first guess that the number of non-zero terms
4712 // will be as many as in b 4749 // will be as many as in b
4713 int x_nz = b.nnz (); 4750 int x_nz = b.nnz ();
4714 int ii = 0; 4751 int ii = 0;
4718 4755
4719 retval.xcidx(0) = 0; 4756 retval.xcidx(0) = 0;
4720 for (int j = 0; j < b_nc; j++) 4757 for (int j = 0; j < b_nc; j++)
4721 { 4758 {
4722 4759
4760 #ifdef UMFPACK_SEPARATE_SPLIT
4723 for (int i = 0; i < b_nr; i++) 4761 for (int i = 0; i < b_nr; i++)
4724 Bx[i] = b.elem (i, j); 4762 Bx[i] = b.elem (i, j);
4725 4763
4726 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, 4764 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai,
4727 X_CAST (const double *, Ax), 4765 X_CAST (const double *, Ax),
4728 NULL, 4766 NULL,
4729 X_CAST (double *, Xx), NULL, 4767 X_CAST (double *, Xx), NULL,
4730 Bx, Bz, Numeric, control, 4768 Bx, Bz, Numeric, control,
4731 info); 4769 info);
4770 #else
4771 for (int i = 0; i < b_nr; i++)
4772 Bz[i] = b.elem (i, j);
4773
4774 status = umfpack_zi_solve (UMFPACK_A, Ap, Ai,
4775 X_CAST (const double *, Ax),
4776 NULL,
4777 X_CAST (double *, Xx), NULL,
4778 X_CAST (double *, Bz), NULL,
4779 Numeric, control,
4780 info);
4781 #endif
4732 if (status < 0) 4782 if (status < 0)
4733 { 4783 {
4734 (*current_liboctave_error_handler) 4784 (*current_liboctave_error_handler)
4735 ("SparseComplexMatrix::solve solve failed"); 4785 ("SparseComplexMatrix::solve solve failed");
4736 4786
4784 4834
4785 umfpack_zi_report_info (control, info); 4835 umfpack_zi_report_info (control, info);
4786 4836
4787 umfpack_zi_free_numeric (&Numeric); 4837 umfpack_zi_free_numeric (&Numeric);
4788 } 4838 }
4839 #else
4840 (*current_liboctave_error_handler) ("UMFPACK not installed");
4841 #endif
4789 } 4842 }
4790 else if (typ != SparseType::Hermitian) 4843 else if (typ != SparseType::Hermitian)
4791 (*current_liboctave_error_handler) ("incorrect matrix type"); 4844 (*current_liboctave_error_handler) ("incorrect matrix type");
4792 } 4845 }
4793 4846
4826 typ = SparseType::Full; 4879 typ = SparseType::Full;
4827 } 4880 }
4828 4881
4829 if (typ == SparseType::Full) 4882 if (typ == SparseType::Full)
4830 { 4883 {
4884 #ifdef HAVE_UMFPACK
4831 Matrix Control, Info; 4885 Matrix Control, Info;
4832 void *Numeric = factorize (err, rcond, Control, Info, sing_handler); 4886 void *Numeric = factorize (err, rcond, Control, Info, sing_handler);
4833 4887
4834 if (err == 0) 4888 if (err == 0)
4835 { 4889 {
4889 4943
4890 umfpack_zi_report_info (control, info); 4944 umfpack_zi_report_info (control, info);
4891 4945
4892 umfpack_zi_free_numeric (&Numeric); 4946 umfpack_zi_free_numeric (&Numeric);
4893 } 4947 }
4948 #else
4949 (*current_liboctave_error_handler) ("UMFPACK not installed");
4950 #endif
4894 } 4951 }
4895 else if (typ != SparseType::Hermitian) 4952 else if (typ != SparseType::Hermitian)
4896 (*current_liboctave_error_handler) ("incorrect matrix type"); 4953 (*current_liboctave_error_handler) ("incorrect matrix type");
4897 } 4954 }
4898 4955
4931 typ = SparseType::Full; 4988 typ = SparseType::Full;
4932 } 4989 }
4933 4990
4934 if (typ == SparseType::Full) 4991 if (typ == SparseType::Full)
4935 { 4992 {
4993 #ifdef HAVE_UMFPACK
4936 Matrix Control, Info; 4994 Matrix Control, Info;
4937 void *Numeric = factorize (err, rcond, Control, Info, sing_handler); 4995 void *Numeric = factorize (err, rcond, Control, Info, sing_handler);
4938 4996
4939 if (err == 0) 4997 if (err == 0)
4940 { 4998 {
5024 5082
5025 umfpack_zi_report_info (control, info); 5083 umfpack_zi_report_info (control, info);
5026 5084
5027 umfpack_zi_free_numeric (&Numeric); 5085 umfpack_zi_free_numeric (&Numeric);
5028 } 5086 }
5087 #else
5088 (*current_liboctave_error_handler) ("UMFPACK not installed");
5089 #endif
5029 } 5090 }
5030 else if (typ != SparseType::Hermitian) 5091 else if (typ != SparseType::Hermitian)
5031 (*current_liboctave_error_handler) ("incorrect matrix type"); 5092 (*current_liboctave_error_handler) ("incorrect matrix type");
5032 } 5093 }
5033 5094