Mercurial > hg > octave-nkf

comparison liboctave/CSparse.cc @ 7482:29980c6b8604

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
don't check f77_exception_encountered
author John W. Eaton <jwe@octave.org>
date Thu, 14 Feb 2008 21:57:50 -0500
parents 402168152bb9
children 8c32f95c2639
comparison
equal deleted inserted replaced
7481:78f3811155f7 7482:29980c6b8604
3747 Complex *result = retval.fortran_vec (); 3747 Complex *result = retval.fortran_vec ();
3748 3748
3749 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result, 3749 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
3750 b.rows(), err)); 3750 b.rows(), err));
3751 3751
3752 if (f77_exception_encountered) 3752 if (err != 0)
3753 (*current_liboctave_error_handler)
3754 ("unrecoverable error in zptsv");
3755 else if (err != 0)
3756 { 3753 {
3757 err = 0; 3754 err = 0;
3758 mattype.mark_as_unsymmetric (); 3755 mattype.mark_as_unsymmetric ();
3759 typ = MatrixType::Tridiagonal; 3756 typ = MatrixType::Tridiagonal;
3760 } 3757 }
3807 Complex *result = retval.fortran_vec (); 3804 Complex *result = retval.fortran_vec ();
3808 3805
3809 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result, 3806 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
3810 b.rows(), err)); 3807 b.rows(), err));
3811 3808
3812 if (f77_exception_encountered) 3809 if (err != 0)
3813 (*current_liboctave_error_handler)
3814 ("unrecoverable error in zgtsv");
3815 else if (err != 0)
3816 { 3810 {
3817 rcond = 0.; 3811 rcond = 0.;
3818 err = -2; 3812 err = -2;
3819 3813
3820 if (sing_handler) 3814 if (sing_handler)
3908 } 3902 }
3909 } 3903 }
3910 3904
3911 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); 3905 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
3912 3906
3913 if (f77_exception_encountered) 3907 if (err != 0)
3914 (*current_liboctave_error_handler) 3908 {
3915 ("unrecoverable error in zgttrf"); 3909 err = -2;
3916 else 3910 rcond = 0.0;
3917 { 3911
3918 if (err != 0) 3912 if (sing_handler)
3919 { 3913 {
3920 err = -2; 3914 sing_handler (rcond);
3921 rcond = 0.0; 3915 mattype.mark_as_rectangular ();
3922 3916 }
3923 if (sing_handler) 3917 else
3924 { 3918 (*current_liboctave_error_handler)
3925 sing_handler (rcond); 3919 ("matrix singular to machine precision");
3926 mattype.mark_as_rectangular (); 3920
3927 } 3921 }
3928 else 3922 else
3929 (*current_liboctave_error_handler) 3923 {
3930 ("matrix singular to machine precision"); 3924 char job = 'N';
3931 3925 volatile octave_idx_type x_nz = b.nnz ();
3932 } 3926 octave_idx_type b_nc = b.cols ();
3933 else 3927 retval = SparseComplexMatrix (nr, b_nc, x_nz);
3934 { 3928 retval.xcidx(0) = 0;
3935 char job = 'N'; 3929 volatile octave_idx_type ii = 0;
3936 volatile octave_idx_type x_nz = b.nnz (); 3930 rcond = 1.0;
3937 octave_idx_type b_nc = b.cols (); 3931
3938 retval = SparseComplexMatrix (nr, b_nc, x_nz); 3932 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
3939 retval.xcidx(0) = 0; 3933
3940 volatile octave_idx_type ii = 0; 3934 for (volatile octave_idx_type j = 0; j < b_nc; j++)
3941 rcond = 1.0; 3935 {
3942 3936 for (octave_idx_type i = 0; i < nr; i++)
3943 OCTAVE_LOCAL_BUFFER (Complex, work, nr); 3937 work[i] = 0.;
3944 3938 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
3945 for (volatile octave_idx_type j = 0; j < b_nc; j++) 3939 work[b.ridx(i)] = b.data(i);
3946 { 3940
3947 for (octave_idx_type i = 0; i < nr; i++) 3941 F77_XFCN (zgttrs, ZGTTRS,
3948 work[i] = 0.; 3942 (F77_CONST_CHAR_ARG2 (&job, 1),
3949 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 3943 nr, 1, DL, D, DU, DU2, pipvt,
3950 work[b.ridx(i)] = b.data(i); 3944 work, b.rows (), err
3951 3945 F77_CHAR_ARG_LEN (1)));
3952 F77_XFCN (zgttrs, ZGTTRS, 3946
3953 (F77_CONST_CHAR_ARG2 (&job, 1), 3947 // Count non-zeros in work vector and adjust
3954 nr, 1, DL, D, DU, DU2, pipvt, 3948 // space in retval if needed
3955 work, b.rows (), err 3949 octave_idx_type new_nnz = 0;
3956 F77_CHAR_ARG_LEN (1))); 3950 for (octave_idx_type i = 0; i < nr; i++)
3957 3951 if (work[i] != 0.)
3958 if (f77_exception_encountered) 3952 new_nnz++;
3959 { 3953
3960 (*current_liboctave_error_handler) 3954 if (ii + new_nnz > x_nz)
3961 ("unrecoverable error in zgttrs"); 3955 {
3962 break; 3956 // Resize the sparse matrix
3963 } 3957 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3964 3958 retval.change_capacity (sz);
3965 // Count non-zeros in work vector and adjust 3959 x_nz = sz;
3966 // space in retval if needed 3960 }
3967 octave_idx_type new_nnz = 0; 3961
3968 for (octave_idx_type i = 0; i < nr; i++) 3962 for (octave_idx_type i = 0; i < nr; i++)
3969 if (work[i] != 0.) 3963 if (work[i] != 0.)
3970 new_nnz++; 3964 {
3971 3965 retval.xridx(ii) = i;
3972 if (ii + new_nnz > x_nz) 3966 retval.xdata(ii++) = work[i];
3973 { 3967 }
3974 // Resize the sparse matrix 3968 retval.xcidx(j+1) = ii;
3975 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; 3969 }
3976 retval.change_capacity (sz); 3970
3977 x_nz = sz; 3971 retval.maybe_compress ();
3978 }
3979
3980 for (octave_idx_type i = 0; i < nr; i++)
3981 if (work[i] != 0.)
3982 {
3983 retval.xridx(ii) = i;
3984 retval.xdata(ii++) = work[i];
3985 }
3986 retval.xcidx(j+1) = ii;
3987 }
3988
3989 retval.maybe_compress ();
3990 }
3991 } 3972 }
3992 } 3973 }
3993 else if (typ != MatrixType::Tridiagonal_Hermitian) 3974 else if (typ != MatrixType::Tridiagonal_Hermitian)
3994 (*current_liboctave_error_handler) ("incorrect matrix type"); 3975 (*current_liboctave_error_handler) ("incorrect matrix type");
3995 } 3976 }
4067 Complex *result = retval.fortran_vec (); 4048 Complex *result = retval.fortran_vec ();
4068 4049
4069 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result, 4050 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
4070 b_nr, err)); 4051 b_nr, err));
4071 4052
4072 if (f77_exception_encountered) 4053 if (err != 0)
4073 {
4074 (*current_liboctave_error_handler)
4075 ("unrecoverable error in zptsv");
4076 err = -1;
4077 }
4078 else if (err != 0)
4079 { 4054 {
4080 err = 0; 4055 err = 0;
4081 mattype.mark_as_unsymmetric (); 4056 mattype.mark_as_unsymmetric ();
4082 typ = MatrixType::Tridiagonal; 4057 typ = MatrixType::Tridiagonal;
4083 } 4058 }
4131 Complex *result = retval.fortran_vec (); 4106 Complex *result = retval.fortran_vec ();
4132 4107
4133 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result, 4108 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
4134 b_nr, err)); 4109 b_nr, err));
4135 4110
4136 if (f77_exception_encountered) 4111 if (err != 0)
4137 {
4138 (*current_liboctave_error_handler)
4139 ("unrecoverable error in zgtsv");
4140 err = -1;
4141 }
4142 else if (err != 0)
4143 { 4112 {
4144 rcond = 0.; 4113 rcond = 0.;
4145 err = -2; 4114 err = -2;
4146 4115
4147 if (sing_handler) 4116 if (sing_handler)
4233 } 4202 }
4234 } 4203 }
4235 4204
4236 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); 4205 F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
4237 4206
4238 if (f77_exception_encountered) 4207 if (err != 0)
4239 (*current_liboctave_error_handler) 4208 {
4240 ("unrecoverable error in zgttrf"); 4209 rcond = 0.0;
4241 else 4210 err = -2;
4242 { 4211
4243 if (err != 0) 4212 if (sing_handler)
4244 { 4213 {
4245 rcond = 0.0; 4214 sing_handler (rcond);
4246 err = -2; 4215 mattype.mark_as_rectangular ();
4247 4216 }
4248 if (sing_handler) 4217 else
4249 { 4218 (*current_liboctave_error_handler)
4250 sing_handler (rcond); 4219 ("matrix singular to machine precision");
4251 mattype.mark_as_rectangular (); 4220 }
4252 } 4221 else
4253 else 4222 {
4254 (*current_liboctave_error_handler) 4223 rcond = 1.;
4255 ("matrix singular to machine precision"); 4224 char job = 'N';
4256 } 4225 octave_idx_type b_nr = b.rows ();
4257 else 4226 octave_idx_type b_nc = b.cols ();
4258 { 4227 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
4259 rcond = 1.; 4228
4260 char job = 'N'; 4229 // Take a first guess that the number of non-zero terms
4261 octave_idx_type b_nr = b.rows (); 4230 // will be as many as in b
4262 octave_idx_type b_nc = b.cols (); 4231 volatile octave_idx_type x_nz = b.nnz ();
4263 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); 4232 volatile octave_idx_type ii = 0;
4264 4233 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4265 // Take a first guess that the number of non-zero terms 4234
4266 // will be as many as in b 4235 retval.xcidx(0) = 0;
4267 volatile octave_idx_type x_nz = b.nnz (); 4236 for (volatile octave_idx_type j = 0; j < b_nc; j++)
4268 volatile octave_idx_type ii = 0; 4237 {
4269 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); 4238
4270 4239 for (octave_idx_type i = 0; i < b_nr; i++)
4271 retval.xcidx(0) = 0; 4240 Bx[i] = b (i,j);
4272 for (volatile octave_idx_type j = 0; j < b_nc; j++) 4241
4273 { 4242 F77_XFCN (zgttrs, ZGTTRS,
4274 4243 (F77_CONST_CHAR_ARG2 (&job, 1),
4275 for (octave_idx_type i = 0; i < b_nr; i++) 4244 nr, 1, DL, D, DU, DU2, pipvt,
4276 Bx[i] = b (i,j); 4245 Bx, b_nr, err
4277 4246 F77_CHAR_ARG_LEN (1)));
4278 F77_XFCN (zgttrs, ZGTTRS, 4247
4279 (F77_CONST_CHAR_ARG2 (&job, 1), 4248 if (err != 0)
4280 nr, 1, DL, D, DU, DU2, pipvt, 4249 {
4281 Bx, b_nr, err 4250 (*current_liboctave_error_handler)
4282 F77_CHAR_ARG_LEN (1))); 4251 ("SparseComplexMatrix::solve solve failed");
4283 4252
4284 if (f77_exception_encountered) 4253 err = -1;
4285 { 4254 break;
4286 (*current_liboctave_error_handler) 4255 }
4287 ("unrecoverable error in zgttrs"); 4256
4288 break; 4257 // Count non-zeros in work vector and adjust
4289 } 4258 // space in retval if needed
4290 4259 octave_idx_type new_nnz = 0;
4291 if (err != 0) 4260 for (octave_idx_type i = 0; i < nr; i++)
4292 { 4261 if (Bx[i] != 0.)
4293 (*current_liboctave_error_handler) 4262 new_nnz++;
4294 ("SparseComplexMatrix::solve solve failed"); 4263
4295 4264 if (ii + new_nnz > x_nz)
4296 err = -1; 4265 {
4297 break; 4266 // Resize the sparse matrix
4298 } 4267 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4299 4268 retval.change_capacity (sz);
4300 // Count non-zeros in work vector and adjust 4269 x_nz = sz;
4301 // space in retval if needed 4270 }
4302 octave_idx_type new_nnz = 0; 4271
4303 for (octave_idx_type i = 0; i < nr; i++) 4272 for (octave_idx_type i = 0; i < nr; i++)
4304 if (Bx[i] != 0.) 4273 if (Bx[i] != 0.)
4305 new_nnz++; 4274 {
4306 4275 retval.xridx(ii) = i;
4307 if (ii + new_nnz > x_nz) 4276 retval.xdata(ii++) = Bx[i];
4308 { 4277 }
4309 // Resize the sparse matrix 4278
4310 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; 4279 retval.xcidx(j+1) = ii;
4311 retval.change_capacity (sz); 4280 }
4312 x_nz = sz; 4281
4313 } 4282 retval.maybe_compress ();
4314
4315 for (octave_idx_type i = 0; i < nr; i++)
4316 if (Bx[i] != 0.)
4317 {
4318 retval.xridx(ii) = i;
4319 retval.xdata(ii++) = Bx[i];
4320 }
4321
4322 retval.xcidx(j+1) = ii;
4323 }
4324
4325 retval.maybe_compress ();
4326 }
4327 } 4283 }
4328 } 4284 }
4329 else if (typ != MatrixType::Tridiagonal_Hermitian) 4285 else if (typ != MatrixType::Tridiagonal_Hermitian)
4330 (*current_liboctave_error_handler) ("incorrect matrix type"); 4286 (*current_liboctave_error_handler) ("incorrect matrix type");
4331 } 4287 }
4388 char job = 'L'; 4344 char job = 'L';
4389 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), 4345 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4390 nr, n_lower, tmp_data, ldm, err 4346 nr, n_lower, tmp_data, ldm, err
4391 F77_CHAR_ARG_LEN (1))); 4347 F77_CHAR_ARG_LEN (1)));
4392 4348
4393 if (f77_exception_encountered) 4349 if (err != 0)
4394 (*current_liboctave_error_handler) 4350 {
4395 ("unrecoverable error in zpbtrf"); 4351 rcond = 0.0;
4396 else 4352 // Matrix is not positive definite!! Fall through to
4397 { 4353 // unsymmetric banded solver.
4398 if (err != 0) 4354 mattype.mark_as_unsymmetric ();
4399 { 4355 typ = MatrixType::Banded;
4400 rcond = 0.0; 4356 err = 0;
4401 // Matrix is not positive definite!! Fall through to 4357 }
4402 // unsymmetric banded solver. 4358 else
4403 mattype.mark_as_unsymmetric (); 4359 {
4404 typ = MatrixType::Banded; 4360 if (calc_cond)
4405 err = 0; 4361 {
4406 } 4362 Array<Complex> z (2 * nr);
4407 else 4363 Complex *pz = z.fortran_vec ();
4408 { 4364 Array<double> iz (nr);
4409 if (calc_cond) 4365 double *piz = iz.fortran_vec ();
4410 { 4366
4411 Array<Complex> z (2 * nr); 4367 F77_XFCN (zpbcon, ZPBCON,
4412 Complex *pz = z.fortran_vec (); 4368 (F77_CONST_CHAR_ARG2 (&job, 1),
4413 Array<double> iz (nr); 4369 nr, n_lower, tmp_data, ldm,
4414 double *piz = iz.fortran_vec (); 4370 anorm, rcond, pz, piz, err
4415 4371 F77_CHAR_ARG_LEN (1)));
4416 F77_XFCN (zpbcon, ZPBCON, 4372
4417 (F77_CONST_CHAR_ARG2 (&job, 1), 4373 if (err != 0)
4418 nr, n_lower, tmp_data, ldm, 4374 err = -2;
4419 anorm, rcond, pz, piz, err 4375
4420 F77_CHAR_ARG_LEN (1))); 4376 volatile double rcond_plus_one = rcond + 1.0;
4421 4377
4422 if (f77_exception_encountered) 4378 if (rcond_plus_one == 1.0 || xisnan (rcond))
4423 (*current_liboctave_error_handler) 4379 {
4424 ("unrecoverable error in zpbcon"); 4380 err = -2;
4425 4381
4426 if (err != 0) 4382 if (sing_handler)
4427 err = -2; 4383 {
4428 4384 sing_handler (rcond);
4429 volatile double rcond_plus_one = rcond + 1.0; 4385 mattype.mark_as_rectangular ();
4430 4386 }
4431 if (rcond_plus_one == 1.0 || xisnan (rcond)) 4387 else
4432 {
4433 err = -2;
4434
4435 if (sing_handler)
4436 {
4437 sing_handler (rcond);
4438 mattype.mark_as_rectangular ();
4439 }
4440 else
4441 (*current_liboctave_error_handler)
4442 ("matrix singular to machine precision, rcond = %g",
4443 rcond);
4444 }
4445 }
4446 else
4447 rcond = 1.0;
4448
4449 if (err == 0)
4450 {
4451 retval = ComplexMatrix (b);
4452 Complex *result = retval.fortran_vec ();
4453
4454 octave_idx_type b_nc = b.cols ();
4455
4456 F77_XFCN (zpbtrs, ZPBTRS,
4457 (F77_CONST_CHAR_ARG2 (&job, 1),
4458 nr, n_lower, b_nc, tmp_data,
4459 ldm, result, b.rows(), err
4460 F77_CHAR_ARG_LEN (1)));
4461
4462 if (f77_exception_encountered)
4463 (*current_liboctave_error_handler) 4388 (*current_liboctave_error_handler)
4464 ("unrecoverable error in zpbtrs"); 4389 ("matrix singular to machine precision, rcond = %g",
4465 4390 rcond);
4466 if (err != 0) 4391 }
4467 { 4392 }
4468 (*current_liboctave_error_handler) 4393 else
4469 ("SparseMatrix::solve solve failed"); 4394 rcond = 1.0;
4470 err = -1; 4395
4471 } 4396 if (err == 0)
4397 {
4398 retval = ComplexMatrix (b);
4399 Complex *result = retval.fortran_vec ();
4400
4401 octave_idx_type b_nc = b.cols ();
4402
4403 F77_XFCN (zpbtrs, ZPBTRS,
4404 (F77_CONST_CHAR_ARG2 (&job, 1),
4405 nr, n_lower, b_nc, tmp_data,
4406 ldm, result, b.rows(), err
4407 F77_CHAR_ARG_LEN (1)));
4408
4409 if (err != 0)
4410 {
4411 (*current_liboctave_error_handler)
4412 ("SparseMatrix::solve solve failed");
4413 err = -1;
4472 } 4414 }
4473 } 4415 }
4474 } 4416 }
4475 } 4417 }
4476 4418
4515 octave_idx_type *pipvt = ipvt.fortran_vec (); 4457 octave_idx_type *pipvt = ipvt.fortran_vec ();
4516 4458
4517 F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper, tmp_data, 4459 F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper, tmp_data,
4518 ldm, pipvt, err)); 4460 ldm, pipvt, err));
4519 4461
4520 if (f77_exception_encountered) 4462 // Throw-away extra info LAPACK gives so as to not
4521 (*current_liboctave_error_handler) 4463 // change output.
4522 ("unrecoverable error in zgbtrf"); 4464 if (err != 0)
4523 else 4465 {
4524 { 4466 rcond = 0.0;
4525 // Throw-away extra info LAPACK gives so as to not 4467 err = -2;
4526 // change output. 4468
4527 if (err != 0) 4469 if (sing_handler)
4528 { 4470 {
4529 rcond = 0.0; 4471 sing_handler (rcond);
4530 err = -2; 4472 mattype.mark_as_rectangular ();
4531 4473 }
4532 if (sing_handler) 4474 else
4533 { 4475 (*current_liboctave_error_handler)
4534 sing_handler (rcond); 4476 ("matrix singular to machine precision");
4535 mattype.mark_as_rectangular (); 4477 }
4536 } 4478 else
4537 else 4479 {
4538 (*current_liboctave_error_handler) 4480 if (calc_cond)
4539 ("matrix singular to machine precision"); 4481 {
4540 } 4482 char job = '1';
4541 else 4483 Array<Complex> z (2 * nr);
4542 { 4484 Complex *pz = z.fortran_vec ();
4543 if (calc_cond) 4485 Array<double> iz (nr);
4544 { 4486 double *piz = iz.fortran_vec ();
4545 char job = '1'; 4487
4546 Array<Complex> z (2 * nr); 4488 F77_XFCN (zgbcon, ZGBCON,
4547 Complex *pz = z.fortran_vec (); 4489 (F77_CONST_CHAR_ARG2 (&job, 1),
4548 Array<double> iz (nr); 4490 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4549 double *piz = iz.fortran_vec (); 4491 anorm, rcond, pz, piz, err
4550 4492 F77_CHAR_ARG_LEN (1)));
4551 F77_XFCN (zgbcon, ZGBCON, 4493
4552 (F77_CONST_CHAR_ARG2 (&job, 1), 4494 if (err != 0)
4553 nc, n_lower, n_upper, tmp_data, ldm, pipvt, 4495 err = -2;
4554 anorm, rcond, pz, piz, err 4496
4555 F77_CHAR_ARG_LEN (1))); 4497 volatile double rcond_plus_one = rcond + 1.0;
4556 4498
4557 if (f77_exception_encountered) 4499 if (rcond_plus_one == 1.0 || xisnan (rcond))
4558 (*current_liboctave_error_handler) 4500 {
4559 ("unrecoverable error in zgbcon"); 4501 err = -2;
4560 4502
4561 if (err != 0) 4503 if (sing_handler)
4562 err = -2; 4504 {
4563 4505 sing_handler (rcond);
4564 volatile double rcond_plus_one = rcond + 1.0; 4506 mattype.mark_as_rectangular ();
4565 4507 }
4566 if (rcond_plus_one == 1.0 || xisnan (rcond)) 4508 else
4567 {
4568 err = -2;
4569
4570 if (sing_handler)
4571 {
4572 sing_handler (rcond);
4573 mattype.mark_as_rectangular ();
4574 }
4575 else
4576 (*current_liboctave_error_handler)
4577 ("matrix singular to machine precision, rcond = %g",
4578 rcond);
4579 }
4580 }
4581 else
4582 rcond = 1.;
4583
4584 if (err == 0)
4585 {
4586 retval = ComplexMatrix (b);
4587 Complex *result = retval.fortran_vec ();
4588
4589 octave_idx_type b_nc = b.cols ();
4590
4591 char job = 'N';
4592 F77_XFCN (zgbtrs, ZGBTRS,
4593 (F77_CONST_CHAR_ARG2 (&job, 1),
4594 nr, n_lower, n_upper, b_nc, tmp_data,
4595 ldm, pipvt, result, b.rows(), err
4596 F77_CHAR_ARG_LEN (1)));
4597
4598 if (f77_exception_encountered)
4599 (*current_liboctave_error_handler) 4509 (*current_liboctave_error_handler)
4600 ("unrecoverable error in zgbtrs"); 4510 ("matrix singular to machine precision, rcond = %g",
4601 } 4511 rcond);
4512 }
4513 }
4514 else
4515 rcond = 1.;
4516
4517 if (err == 0)
4518 {
4519 retval = ComplexMatrix (b);
4520 Complex *result = retval.fortran_vec ();
4521
4522 octave_idx_type b_nc = b.cols ();
4523
4524 char job = 'N';
4525 F77_XFCN (zgbtrs, ZGBTRS,
4526 (F77_CONST_CHAR_ARG2 (&job, 1),
4527 nr, n_lower, n_upper, b_nc, tmp_data,
4528 ldm, pipvt, result, b.rows(), err
4529 F77_CHAR_ARG_LEN (1)));
4602 } 4530 }
4603 } 4531 }
4604 } 4532 }
4605 else if (typ != MatrixType::Banded_Hermitian) 4533 else if (typ != MatrixType::Banded_Hermitian)
4606 (*current_liboctave_error_handler) ("incorrect matrix type"); 4534 (*current_liboctave_error_handler) ("incorrect matrix type");
4665 char job = 'L'; 4593 char job = 'L';
4666 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), 4594 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4667 nr, n_lower, tmp_data, ldm, err 4595 nr, n_lower, tmp_data, ldm, err
4668 F77_CHAR_ARG_LEN (1))); 4596 F77_CHAR_ARG_LEN (1)));
4669 4597
4670 if (f77_exception_encountered) 4598 if (err != 0)
4671 (*current_liboctave_error_handler) 4599 {
4672 ("unrecoverable error in zpbtrf"); 4600 rcond = 0.0;
4673 else 4601 mattype.mark_as_unsymmetric ();
4674 { 4602 typ = MatrixType::Banded;
4675 if (err != 0) 4603 err = 0;
4676 { 4604 }
4677 rcond = 0.0; 4605 else
4678 mattype.mark_as_unsymmetric (); 4606 {
4679 typ = MatrixType::Banded; 4607 if (calc_cond)
4680 err = 0; 4608 {
4681 } 4609 Array<Complex> z (2 * nr);
4682 else 4610 Complex *pz = z.fortran_vec ();
4683 { 4611 Array<double> iz (nr);
4684 if (calc_cond) 4612 double *piz = iz.fortran_vec ();
4685 { 4613
4686 Array<Complex> z (2 * nr); 4614 F77_XFCN (zpbcon, ZPBCON,
4687 Complex *pz = z.fortran_vec (); 4615 (F77_CONST_CHAR_ARG2 (&job, 1),
4688 Array<double> iz (nr); 4616 nr, n_lower, tmp_data, ldm,
4689 double *piz = iz.fortran_vec (); 4617 anorm, rcond, pz, piz, err
4690 4618 F77_CHAR_ARG_LEN (1)));
4691 F77_XFCN (zpbcon, ZPBCON, 4619
4692 (F77_CONST_CHAR_ARG2 (&job, 1), 4620 if (err != 0)
4693 nr, n_lower, tmp_data, ldm, 4621 err = -2;
4694 anorm, rcond, pz, piz, err 4622
4695 F77_CHAR_ARG_LEN (1))); 4623 volatile double rcond_plus_one = rcond + 1.0;
4696 4624
4697 if (f77_exception_encountered) 4625 if (rcond_plus_one == 1.0 || xisnan (rcond))
4698 (*current_liboctave_error_handler) 4626 {
4699 ("unrecoverable error in zpbcon"); 4627 err = -2;
4700 4628
4701 if (err != 0) 4629 if (sing_handler)
4702 err = -2; 4630 {
4703 4631 sing_handler (rcond);
4704 volatile double rcond_plus_one = rcond + 1.0; 4632 mattype.mark_as_rectangular ();
4705 4633 }
4706 if (rcond_plus_one == 1.0 || xisnan (rcond)) 4634 else
4707 { 4635 (*current_liboctave_error_handler)
4708 err = -2; 4636 ("matrix singular to machine precision, rcond = %g",
4709 4637 rcond);
4710 if (sing_handler) 4638 }
4639 }
4640 else
4641 rcond = 1.0;
4642
4643 if (err == 0)
4644 {
4645 octave_idx_type b_nr = b.rows ();
4646 octave_idx_type b_nc = b.cols ();
4647 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
4648
4649 // Take a first guess that the number of non-zero terms
4650 // will be as many as in b
4651 volatile octave_idx_type x_nz = b.nnz ();
4652 volatile octave_idx_type ii = 0;
4653 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4654
4655 retval.xcidx(0) = 0;
4656 for (volatile octave_idx_type j = 0; j < b_nc; j++)
4657 {
4658 for (octave_idx_type i = 0; i < b_nr; i++)
4659 Bx[i] = b.elem (i, j);
4660
4661 F77_XFCN (zpbtrs, ZPBTRS,
4662 (F77_CONST_CHAR_ARG2 (&job, 1),
4663 nr, n_lower, 1, tmp_data,
4664 ldm, Bx, b_nr, err
4665 F77_CHAR_ARG_LEN (1)));
4666
4667 if (err != 0)
4668 {
4669 (*current_liboctave_error_handler)
4670 ("SparseComplexMatrix::solve solve failed");
4671 err = -1;
4672 break;
4673 }
4674
4675 for (octave_idx_type i = 0; i < b_nr; i++)
4676 {
4677 Complex tmp = Bx[i];
4678 if (tmp != 0.0)
4711 { 4679 {
4712 sing_handler (rcond); 4680 if (ii == x_nz)
4713 mattype.mark_as_rectangular (); 4681 {
4682 // Resize the sparse matrix
4683 octave_idx_type sz = x_nz *
4684 (b_nc - j) / b_nc;
4685 sz = (sz > 10 ? sz : 10) + x_nz;
4686 retval.change_capacity (sz);
4687 x_nz = sz;
4688 }
4689 retval.xdata(ii) = tmp;
4690 retval.xridx(ii++) = i;
4714 } 4691 }
4715 else
4716 (*current_liboctave_error_handler)
4717 ("matrix singular to machine precision, rcond = %g",
4718 rcond);
4719 }
4720 }
4721 else
4722 rcond = 1.0;
4723
4724 if (err == 0)
4725 {
4726 octave_idx_type b_nr = b.rows ();
4727 octave_idx_type b_nc = b.cols ();
4728 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
4729
4730 // Take a first guess that the number of non-zero terms
4731 // will be as many as in b
4732 volatile octave_idx_type x_nz = b.nnz ();
4733 volatile octave_idx_type ii = 0;
4734 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4735
4736 retval.xcidx(0) = 0;
4737 for (volatile octave_idx_type j = 0; j < b_nc; j++)
4738 {
4739 for (octave_idx_type i = 0; i < b_nr; i++)
4740 Bx[i] = b.elem (i, j);
4741
4742 F77_XFCN (zpbtrs, ZPBTRS,
4743 (F77_CONST_CHAR_ARG2 (&job, 1),
4744 nr, n_lower, 1, tmp_data,
4745 ldm, Bx, b_nr, err
4746 F77_CHAR_ARG_LEN (1)));
4747
4748 if (f77_exception_encountered)
4749 {
4750 (*current_liboctave_error_handler)
4751 ("unrecoverable error in dpbtrs");
4752 err = -1;
4753 break;
4754 }
4755
4756 if (err != 0)
4757 {
4758 (*current_liboctave_error_handler)
4759 ("SparseComplexMatrix::solve solve failed");
4760 err = -1;
4761 break;
4762 }
4763
4764 for (octave_idx_type i = 0; i < b_nr; i++)
4765 {
4766 Complex tmp = Bx[i];
4767 if (tmp != 0.0)
4768 {
4769 if (ii == x_nz)
4770 {
4771 // Resize the sparse matrix
4772 octave_idx_type sz = x_nz *
4773 (b_nc - j) / b_nc;
4774 sz = (sz > 10 ? sz : 10) + x_nz;
4775 retval.change_capacity (sz);
4776 x_nz = sz;
4777 }
4778 retval.xdata(ii) = tmp;
4779 retval.xridx(ii++) = i;
4780 }
4781 }
4782 retval.xcidx(j+1) = ii;
4783 } 4692 }
4784 4693 retval.xcidx(j+1) = ii;
4785 retval.maybe_compress (); 4694 }
4786 } 4695
4696 retval.maybe_compress ();
4787 } 4697 }
4788 } 4698 }
4789 } 4699 }
4790 4700
4791 if (typ == MatrixType::Banded) 4701 if (typ == MatrixType::Banded)
4829 octave_idx_type *pipvt = ipvt.fortran_vec (); 4739 octave_idx_type *pipvt = ipvt.fortran_vec ();
4830 4740
4831 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, 4741 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
4832 ldm, pipvt, err)); 4742 ldm, pipvt, err));
4833 4743
4834 if (f77_exception_encountered) 4744 if (err != 0)
4835 (*current_liboctave_error_handler) 4745 {
4836 ("unrecoverable error in zgbtrf"); 4746 rcond = 0.0;
4837 else 4747 err = -2;
4838 { 4748
4839 if (err != 0) 4749 if (sing_handler)
4840 { 4750 {
4841 rcond = 0.0; 4751 sing_handler (rcond);
4842 err = -2; 4752 mattype.mark_as_rectangular ();
4843 4753 }
4844 if (sing_handler) 4754 else
4845 { 4755 (*current_liboctave_error_handler)
4846 sing_handler (rcond); 4756 ("matrix singular to machine precision");
4847 mattype.mark_as_rectangular (); 4757
4848 } 4758 }
4849 else 4759 else
4850 (*current_liboctave_error_handler) 4760 {
4851 ("matrix singular to machine precision"); 4761 if (calc_cond)
4852 4762 {
4853 } 4763 char job = '1';
4854 else 4764 Array<Complex> z (2 * nr);
4855 { 4765 Complex *pz = z.fortran_vec ();
4856 if (calc_cond) 4766 Array<double> iz (nr);
4857 { 4767 double *piz = iz.fortran_vec ();
4858 char job = '1'; 4768
4859 Array<Complex> z (2 * nr); 4769 F77_XFCN (zgbcon, ZGBCON,
4860 Complex *pz = z.fortran_vec (); 4770 (F77_CONST_CHAR_ARG2 (&job, 1),
4861 Array<double> iz (nr); 4771 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4862 double *piz = iz.fortran_vec (); 4772 anorm, rcond, pz, piz, err
4863 4773 F77_CHAR_ARG_LEN (1)));
4864 F77_XFCN (zgbcon, ZGBCON, 4774
4865 (F77_CONST_CHAR_ARG2 (&job, 1), 4775 if (err != 0)
4866 nc, n_lower, n_upper, tmp_data, ldm, pipvt, 4776 err = -2;
4867 anorm, rcond, pz, piz, err 4777
4868 F77_CHAR_ARG_LEN (1))); 4778 volatile double rcond_plus_one = rcond + 1.0;
4869 4779
4870 if (f77_exception_encountered) 4780 if (rcond_plus_one == 1.0 || xisnan (rcond))
4871 (*current_liboctave_error_handler) 4781 {
4872 ("unrecoverable error in zgbcon"); 4782 err = -2;
4873 4783
4874 if (err != 0) 4784 if (sing_handler)
4875 err = -2;
4876
4877 volatile double rcond_plus_one = rcond + 1.0;
4878
4879 if (rcond_plus_one == 1.0 || xisnan (rcond))
4880 {
4881 err = -2;
4882
4883 if (sing_handler)
4884 {
4885 sing_handler (rcond);
4886 mattype.mark_as_rectangular ();
4887 }
4888 else
4889 (*current_liboctave_error_handler)
4890 ("matrix singular to machine precision, rcond = %g",
4891 rcond);
4892 }
4893 }
4894 else
4895 rcond = 1.;
4896
4897 if (err == 0)
4898 {
4899 char job = 'N';
4900 volatile octave_idx_type x_nz = b.nnz ();
4901 octave_idx_type b_nc = b.cols ();
4902 retval = SparseComplexMatrix (nr, b_nc, x_nz);
4903 retval.xcidx(0) = 0;
4904 volatile octave_idx_type ii = 0;
4905
4906 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
4907
4908 for (volatile octave_idx_type j = 0; j < b_nc; j++)
4909 { 4785 {
4910 for (octave_idx_type i = 0; i < nr; i++) 4786 sing_handler (rcond);
4911 work[i] = 0.; 4787 mattype.mark_as_rectangular ();
4912 for (octave_idx_type i = b.cidx(j);
4913 i < b.cidx(j+1); i++)
4914 work[b.ridx(i)] = b.data(i);
4915
4916 F77_XFCN (zgbtrs, ZGBTRS,
4917 (F77_CONST_CHAR_ARG2 (&job, 1),
4918 nr, n_lower, n_upper, 1, tmp_data,
4919 ldm, pipvt, work, b.rows (), err
4920 F77_CHAR_ARG_LEN (1)));
4921
4922 if (f77_exception_encountered)
4923 {
4924 (*current_liboctave_error_handler)
4925 ("unrecoverable error in zgbtrs");
4926 break;
4927 }
4928
4929 // Count non-zeros in work vector and adjust
4930 // space in retval if needed
4931 octave_idx_type new_nnz = 0;
4932 for (octave_idx_type i = 0; i < nr; i++)
4933 if (work[i] != 0.)
4934 new_nnz++;
4935
4936 if (ii + new_nnz > x_nz)
4937 {
4938 // Resize the sparse matrix
4939 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4940 retval.change_capacity (sz);
4941 x_nz = sz;
4942 }
4943
4944 for (octave_idx_type i = 0; i < nr; i++)
4945 if (work[i] != 0.)
4946 {
4947 retval.xridx(ii) = i;
4948 retval.xdata(ii++) = work[i];
4949 }
4950 retval.xcidx(j+1) = ii;
4951 } 4788 }
4952 4789 else
4953 retval.maybe_compress (); 4790 (*current_liboctave_error_handler)
4954 } 4791 ("matrix singular to machine precision, rcond = %g",
4792 rcond);
4793 }
4794 }
4795 else
4796 rcond = 1.;
4797
4798 if (err == 0)
4799 {
4800 char job = 'N';
4801 volatile octave_idx_type x_nz = b.nnz ();
4802 octave_idx_type b_nc = b.cols ();
4803 retval = SparseComplexMatrix (nr, b_nc, x_nz);
4804 retval.xcidx(0) = 0;
4805 volatile octave_idx_type ii = 0;
4806
4807 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
4808
4809 for (volatile octave_idx_type j = 0; j < b_nc; j++)
4810 {
4811 for (octave_idx_type i = 0; i < nr; i++)
4812 work[i] = 0.;
4813 for (octave_idx_type i = b.cidx(j);
4814 i < b.cidx(j+1); i++)
4815 work[b.ridx(i)] = b.data(i);
4816
4817 F77_XFCN (zgbtrs, ZGBTRS,
4818 (F77_CONST_CHAR_ARG2 (&job, 1),
4819 nr, n_lower, n_upper, 1, tmp_data,
4820 ldm, pipvt, work, b.rows (), err
4821 F77_CHAR_ARG_LEN (1)));
4822
4823 // Count non-zeros in work vector and adjust
4824 // space in retval if needed
4825 octave_idx_type new_nnz = 0;
4826 for (octave_idx_type i = 0; i < nr; i++)
4827 if (work[i] != 0.)
4828 new_nnz++;
4829
4830 if (ii + new_nnz > x_nz)
4831 {
4832 // Resize the sparse matrix
4833 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4834 retval.change_capacity (sz);
4835 x_nz = sz;
4836 }
4837
4838 for (octave_idx_type i = 0; i < nr; i++)
4839 if (work[i] != 0.)
4840 {
4841 retval.xridx(ii) = i;
4842 retval.xdata(ii++) = work[i];
4843 }
4844 retval.xcidx(j+1) = ii;
4845 }
4846
4847 retval.maybe_compress ();
4955 } 4848 }
4956 } 4849 }
4957 } 4850 }
4958 else if (typ != MatrixType::Banded_Hermitian) 4851 else if (typ != MatrixType::Banded_Hermitian)
4959 (*current_liboctave_error_handler) ("incorrect matrix type"); 4852 (*current_liboctave_error_handler) ("incorrect matrix type");
5018 char job = 'L'; 4911 char job = 'L';
5019 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), 4912 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5020 nr, n_lower, tmp_data, ldm, err 4913 nr, n_lower, tmp_data, ldm, err
5021 F77_CHAR_ARG_LEN (1))); 4914 F77_CHAR_ARG_LEN (1)));
5022 4915
5023 if (f77_exception_encountered) 4916 if (err != 0)
5024 (*current_liboctave_error_handler) 4917 {
5025 ("unrecoverable error in zpbtrf"); 4918 // Matrix is not positive definite!! Fall through to
5026 else 4919 // unsymmetric banded solver.
5027 { 4920 rcond = 0.0;
5028 if (err != 0) 4921 mattype.mark_as_unsymmetric ();
5029 { 4922 typ = MatrixType::Banded;
5030 // Matrix is not positive definite!! Fall through to 4923 err = 0;
5031 // unsymmetric banded solver. 4924 }
5032 rcond = 0.0; 4925 else
5033 mattype.mark_as_unsymmetric (); 4926 {
5034 typ = MatrixType::Banded; 4927 if (calc_cond)
5035 err = 0; 4928 {
5036 } 4929 Array<Complex> z (2 * nr);
5037 else 4930 Complex *pz = z.fortran_vec ();
5038 { 4931 Array<double> iz (nr);
5039 if (calc_cond) 4932 double *piz = iz.fortran_vec ();
5040 { 4933
5041 Array<Complex> z (2 * nr); 4934 F77_XFCN (zpbcon, ZPBCON,
5042 Complex *pz = z.fortran_vec (); 4935 (F77_CONST_CHAR_ARG2 (&job, 1),
5043 Array<double> iz (nr); 4936 nr, n_lower, tmp_data, ldm,
5044 double *piz = iz.fortran_vec (); 4937 anorm, rcond, pz, piz, err
5045 4938 F77_CHAR_ARG_LEN (1)));
5046 F77_XFCN (zpbcon, ZPBCON, 4939
5047 (F77_CONST_CHAR_ARG2 (&job, 1), 4940 if (err != 0)
5048 nr, n_lower, tmp_data, ldm, 4941 err = -2;
5049 anorm, rcond, pz, piz, err 4942
5050 F77_CHAR_ARG_LEN (1))); 4943 volatile double rcond_plus_one = rcond + 1.0;
5051 4944
5052 if (f77_exception_encountered) 4945 if (rcond_plus_one == 1.0 || xisnan (rcond))
5053 (*current_liboctave_error_handler) 4946 {
5054 ("unrecoverable error in zpbcon"); 4947 err = -2;
5055 4948
5056 if (err != 0) 4949 if (sing_handler)
5057 err = -2;
5058
5059 volatile double rcond_plus_one = rcond + 1.0;
5060
5061 if (rcond_plus_one == 1.0 || xisnan (rcond))
5062 {
5063 err = -2;
5064
5065 if (sing_handler)
5066 {
5067 sing_handler (rcond);
5068 mattype.mark_as_rectangular ();
5069 }
5070 else
5071 (*current_liboctave_error_handler)
5072 ("matrix singular to machine precision, rcond = %g",
5073 rcond);
5074 }
5075 }
5076 else
5077 rcond = 1.0;
5078
5079 if (err == 0)
5080 {
5081 octave_idx_type b_nr = b.rows ();
5082 octave_idx_type b_nc = b.cols ();
5083 retval = ComplexMatrix (b);
5084 Complex *result = retval.fortran_vec ();
5085
5086 F77_XFCN (zpbtrs, ZPBTRS,
5087 (F77_CONST_CHAR_ARG2 (&job, 1),
5088 nr, n_lower, b_nc, tmp_data,
5089 ldm, result, b_nr, err
5090 F77_CHAR_ARG_LEN (1)));
5091
5092 if (f77_exception_encountered)
5093 { 4950 {
5094 (*current_liboctave_error_handler) 4951 sing_handler (rcond);
5095 ("unrecoverable error in zpbtrs"); 4952 mattype.mark_as_rectangular ();
5096 err = -1;
5097 } 4953 }
5098 4954 else
5099 if (err != 0) 4955 (*current_liboctave_error_handler)
5100 { 4956 ("matrix singular to machine precision, rcond = %g",
5101 (*current_liboctave_error_handler) 4957 rcond);
5102 ("SparseComplexMatrix::solve solve failed"); 4958 }
5103 err = -1; 4959 }
5104 } 4960 else
4961 rcond = 1.0;
4962
4963 if (err == 0)
4964 {
4965 octave_idx_type b_nr = b.rows ();
4966 octave_idx_type b_nc = b.cols ();
4967 retval = ComplexMatrix (b);
4968 Complex *result = retval.fortran_vec ();
4969
4970 F77_XFCN (zpbtrs, ZPBTRS,
4971 (F77_CONST_CHAR_ARG2 (&job, 1),
4972 nr, n_lower, b_nc, tmp_data,
4973 ldm, result, b_nr, err
4974 F77_CHAR_ARG_LEN (1)));
4975
4976 if (err != 0)
4977 {
4978 (*current_liboctave_error_handler)
4979 ("SparseComplexMatrix::solve solve failed");
4980 err = -1;
5105 } 4981 }
5106 } 4982 }
5107 } 4983 }
5108 } 4984 }
5109 4985
5148 octave_idx_type *pipvt = ipvt.fortran_vec (); 5024 octave_idx_type *pipvt = ipvt.fortran_vec ();
5149 5025
5150 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, 5026 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
5151 ldm, pipvt, err)); 5027 ldm, pipvt, err));
5152 5028
5153 if (f77_exception_encountered) 5029 if (err != 0)
5154 (*current_liboctave_error_handler) 5030 {
5155 ("unrecoverable error in zgbtrf"); 5031 err = -2;
5156 else 5032 rcond = 0.0;
5157 { 5033
5158 if (err != 0) 5034 if (sing_handler)
5159 { 5035 {
5160 err = -2; 5036 sing_handler (rcond);
5161 rcond = 0.0; 5037 mattype.mark_as_rectangular ();
5162 5038 }
5163 if (sing_handler) 5039 else
5164 { 5040 (*current_liboctave_error_handler)
5165 sing_handler (rcond); 5041 ("matrix singular to machine precision");
5166 mattype.mark_as_rectangular (); 5042 }
5167 } 5043 else
5168 else 5044 {
5169 (*current_liboctave_error_handler) 5045 if (calc_cond)
5170 ("matrix singular to machine precision"); 5046 {
5171 } 5047 char job = '1';
5172 else 5048 Array<Complex> z (2 * nr);
5173 { 5049 Complex *pz = z.fortran_vec ();
5174 if (calc_cond) 5050 Array<double> iz (nr);
5175 { 5051 double *piz = iz.fortran_vec ();
5176 char job = '1'; 5052
5177 Array<Complex> z (2 * nr); 5053 F77_XFCN (zgbcon, ZGBCON,
5178 Complex *pz = z.fortran_vec (); 5054 (F77_CONST_CHAR_ARG2 (&job, 1),
5179 Array<double> iz (nr); 5055 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5180 double *piz = iz.fortran_vec (); 5056 anorm, rcond, pz, piz, err
5181 5057 F77_CHAR_ARG_LEN (1)));
5182 F77_XFCN (zgbcon, ZGBCON, 5058
5183 (F77_CONST_CHAR_ARG2 (&job, 1), 5059 if (err != 0)
5184 nc, n_lower, n_upper, tmp_data, ldm, pipvt, 5060 err = -2;
5185 anorm, rcond, pz, piz, err 5061
5186 F77_CHAR_ARG_LEN (1))); 5062 volatile double rcond_plus_one = rcond + 1.0;
5187 5063
5188 if (f77_exception_encountered) 5064 if (rcond_plus_one == 1.0 || xisnan (rcond))
5189 (*current_liboctave_error_handler) 5065 {
5190 ("unrecoverable error in zgbcon"); 5066 err = -2;
5191 5067
5192 if (err != 0) 5068 if (sing_handler)
5193 err = -2;
5194
5195 volatile double rcond_plus_one = rcond + 1.0;
5196
5197 if (rcond_plus_one == 1.0 || xisnan (rcond))
5198 {
5199 err = -2;
5200
5201 if (sing_handler)
5202 {
5203 sing_handler (rcond);
5204 mattype.mark_as_rectangular ();
5205 }
5206 else
5207 (*current_liboctave_error_handler)
5208 ("matrix singular to machine precision, rcond = %g",
5209 rcond);
5210 }
5211 }
5212 else
5213 rcond = 1.;
5214
5215 if (err == 0)
5216 {
5217 char job = 'N';
5218 octave_idx_type b_nc = b.cols ();
5219 retval = ComplexMatrix (b);
5220 Complex *result = retval.fortran_vec ();
5221
5222 F77_XFCN (zgbtrs, ZGBTRS,
5223 (F77_CONST_CHAR_ARG2 (&job, 1),
5224 nr, n_lower, n_upper, b_nc, tmp_data,
5225 ldm, pipvt, result, b.rows (), err
5226 F77_CHAR_ARG_LEN (1)));
5227
5228 if (f77_exception_encountered)
5229 { 5069 {
5230 (*current_liboctave_error_handler) 5070 sing_handler (rcond);
5231 ("unrecoverable error in dgbtrs"); 5071 mattype.mark_as_rectangular ();
5232 } 5072 }
5233 } 5073 else
5074 (*current_liboctave_error_handler)
5075 ("matrix singular to machine precision, rcond = %g",
5076 rcond);
5077 }
5078 }
5079 else
5080 rcond = 1.;
5081
5082 if (err == 0)
5083 {
5084 char job = 'N';
5085 octave_idx_type b_nc = b.cols ();
5086 retval = ComplexMatrix (b);
5087 Complex *result = retval.fortran_vec ();
5088
5089 F77_XFCN (zgbtrs, ZGBTRS,
5090 (F77_CONST_CHAR_ARG2 (&job, 1),
5091 nr, n_lower, n_upper, b_nc, tmp_data,
5092 ldm, pipvt, result, b.rows (), err
5093 F77_CHAR_ARG_LEN (1)));
5234 } 5094 }
5235 } 5095 }
5236 } 5096 }
5237 else if (typ != MatrixType::Banded_Hermitian) 5097 else if (typ != MatrixType::Banded_Hermitian)
5238 (*current_liboctave_error_handler) ("incorrect matrix type"); 5098 (*current_liboctave_error_handler) ("incorrect matrix type");
5297 char job = 'L'; 5157 char job = 'L';
5298 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), 5158 F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5299 nr, n_lower, tmp_data, ldm, err 5159 nr, n_lower, tmp_data, ldm, err
5300 F77_CHAR_ARG_LEN (1))); 5160 F77_CHAR_ARG_LEN (1)));
5301 5161
5302 if (f77_exception_encountered) 5162 if (err != 0)
5303 (*current_liboctave_error_handler) 5163 {
5304 ("unrecoverable error in zpbtrf"); 5164 // Matrix is not positive definite!! Fall through to
5305 else 5165 // unsymmetric banded solver.
5306 { 5166 mattype.mark_as_unsymmetric ();
5307 if (err != 0) 5167 typ = MatrixType::Banded;
5308 { 5168
5309 // Matrix is not positive definite!! Fall through to 5169 rcond = 0.0;
5310 // unsymmetric banded solver. 5170 err = 0;
5311 mattype.mark_as_unsymmetric (); 5171 }
5312 typ = MatrixType::Banded; 5172 else
5313 5173 {
5314 rcond = 0.0; 5174 if (calc_cond)
5315 err = 0; 5175 {
5316 } 5176 Array<Complex> z (2 * nr);
5317 else 5177 Complex *pz = z.fortran_vec ();
5318 { 5178 Array<double> iz (nr);
5319 if (calc_cond) 5179 double *piz = iz.fortran_vec ();
5320 { 5180
5321 Array<Complex> z (2 * nr); 5181 F77_XFCN (zpbcon, ZPBCON,
5322 Complex *pz = z.fortran_vec (); 5182 (F77_CONST_CHAR_ARG2 (&job, 1),
5323 Array<double> iz (nr); 5183 nr, n_lower, tmp_data, ldm,
5324 double *piz = iz.fortran_vec (); 5184 anorm, rcond, pz, piz, err
5325 5185 F77_CHAR_ARG_LEN (1)));
5326 F77_XFCN (zpbcon, ZPBCON, 5186
5327 (F77_CONST_CHAR_ARG2 (&job, 1), 5187 if (err != 0)
5328 nr, n_lower, tmp_data, ldm, 5188 err = -2;
5329 anorm, rcond, pz, piz, err 5189
5330 F77_CHAR_ARG_LEN (1))); 5190 volatile double rcond_plus_one = rcond + 1.0;
5331 5191
5332 if (f77_exception_encountered) 5192 if (rcond_plus_one == 1.0 || xisnan (rcond))
5333 (*current_liboctave_error_handler) 5193 {
5334 ("unrecoverable error in zpbcon"); 5194 err = -2;
5335 5195
5336 if (err != 0) 5196 if (sing_handler)
5337 err = -2;
5338
5339 volatile double rcond_plus_one = rcond + 1.0;
5340
5341 if (rcond_plus_one == 1.0 || xisnan (rcond))
5342 {
5343 err = -2;
5344
5345 if (sing_handler)
5346 {
5347 sing_handler (rcond);
5348 mattype.mark_as_rectangular ();
5349 }
5350 else
5351 (*current_liboctave_error_handler)
5352 ("matrix singular to machine precision, rcond = %g",
5353 rcond);
5354 }
5355 }
5356 else
5357 rcond = 1.0;
5358
5359 if (err == 0)
5360 {
5361 octave_idx_type b_nr = b.rows ();
5362 octave_idx_type b_nc = b.cols ();
5363 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
5364
5365 // Take a first guess that the number of non-zero terms
5366 // will be as many as in b
5367 volatile octave_idx_type x_nz = b.nnz ();
5368 volatile octave_idx_type ii = 0;
5369 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
5370
5371 retval.xcidx(0) = 0;
5372 for (volatile octave_idx_type j = 0; j < b_nc; j++)
5373 { 5197 {
5374 5198 sing_handler (rcond);
5375 for (octave_idx_type i = 0; i < b_nr; i++) 5199 mattype.mark_as_rectangular ();
5376 Bx[i] = b (i,j);
5377
5378 F77_XFCN (zpbtrs, ZPBTRS,
5379 (F77_CONST_CHAR_ARG2 (&job, 1),
5380 nr, n_lower, 1, tmp_data,
5381 ldm, Bx, b_nr, err
5382 F77_CHAR_ARG_LEN (1)));
5383
5384 if (f77_exception_encountered)
5385 {
5386 (*current_liboctave_error_handler)
5387 ("unrecoverable error in zpbtrs");
5388 err = -1;
5389 break;
5390 }
5391
5392 if (err != 0)
5393 {
5394 (*current_liboctave_error_handler)
5395 ("SparseMatrix::solve solve failed");
5396 err = -1;
5397 break;
5398 }
5399
5400 // Count non-zeros in work vector and adjust
5401 // space in retval if needed
5402 octave_idx_type new_nnz = 0;
5403 for (octave_idx_type i = 0; i < nr; i++)
5404 if (Bx[i] != 0.)
5405 new_nnz++;
5406
5407 if (ii + new_nnz > x_nz)
5408 {
5409 // Resize the sparse matrix
5410 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5411 retval.change_capacity (sz);
5412 x_nz = sz;
5413 }
5414
5415 for (octave_idx_type i = 0; i < nr; i++)
5416 if (Bx[i] != 0.)
5417 {
5418 retval.xridx(ii) = i;
5419 retval.xdata(ii++) = Bx[i];
5420 }
5421
5422 retval.xcidx(j+1) = ii;
5423 } 5200 }
5424 5201 else
5425 retval.maybe_compress (); 5202 (*current_liboctave_error_handler)
5426 } 5203 ("matrix singular to machine precision, rcond = %g",
5204 rcond);
5205 }
5206 }
5207 else
5208 rcond = 1.0;
5209
5210 if (err == 0)
5211 {
5212 octave_idx_type b_nr = b.rows ();
5213 octave_idx_type b_nc = b.cols ();
5214 OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
5215
5216 // Take a first guess that the number of non-zero terms
5217 // will be as many as in b
5218 volatile octave_idx_type x_nz = b.nnz ();
5219 volatile octave_idx_type ii = 0;
5220 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
5221
5222 retval.xcidx(0) = 0;
5223 for (volatile octave_idx_type j = 0; j < b_nc; j++)
5224 {
5225
5226 for (octave_idx_type i = 0; i < b_nr; i++)
5227 Bx[i] = b (i,j);
5228
5229 F77_XFCN (zpbtrs, ZPBTRS,
5230 (F77_CONST_CHAR_ARG2 (&job, 1),
5231 nr, n_lower, 1, tmp_data,
5232 ldm, Bx, b_nr, err
5233 F77_CHAR_ARG_LEN (1)));
5234
5235 if (err != 0)
5236 {
5237 (*current_liboctave_error_handler)
5238 ("SparseMatrix::solve solve failed");
5239 err = -1;
5240 break;
5241 }
5242
5243 // Count non-zeros in work vector and adjust
5244 // space in retval if needed
5245 octave_idx_type new_nnz = 0;
5246 for (octave_idx_type i = 0; i < nr; i++)
5247 if (Bx[i] != 0.)
5248 new_nnz++;
5249
5250 if (ii + new_nnz > x_nz)
5251 {
5252 // Resize the sparse matrix
5253 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5254 retval.change_capacity (sz);
5255 x_nz = sz;
5256 }
5257
5258 for (octave_idx_type i = 0; i < nr; i++)
5259 if (Bx[i] != 0.)
5260 {
5261 retval.xridx(ii) = i;
5262 retval.xdata(ii++) = Bx[i];
5263 }
5264
5265 retval.xcidx(j+1) = ii;
5266 }
5267
5268 retval.maybe_compress ();
5427 } 5269 }
5428 } 5270 }
5429 } 5271 }
5430 5272
5431 if (typ == MatrixType::Banded) 5273 if (typ == MatrixType::Banded)
5469 octave_idx_type *pipvt = ipvt.fortran_vec (); 5311 octave_idx_type *pipvt = ipvt.fortran_vec ();
5470 5312
5471 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, 5313 F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
5472 ldm, pipvt, err)); 5314 ldm, pipvt, err));
5473 5315
5474 if (f77_exception_encountered) 5316 if (err != 0)
5475 (*current_liboctave_error_handler) 5317 {
5476 ("unrecoverable error in xgbtrf"); 5318 err = -2;
5477 else 5319 rcond = 0.0;
5478 { 5320
5479 if (err != 0) 5321 if (sing_handler)
5480 { 5322 {
5481 err = -2; 5323 sing_handler (rcond);
5482 rcond = 0.0; 5324 mattype.mark_as_rectangular ();
5483 5325 }
5484 if (sing_handler) 5326 else
5485 { 5327 (*current_liboctave_error_handler)
5486 sing_handler (rcond); 5328 ("matrix singular to machine precision");
5487 mattype.mark_as_rectangular (); 5329
5488 } 5330 }
5489 else 5331 else
5490 (*current_liboctave_error_handler) 5332 {
5491 ("matrix singular to machine precision"); 5333 if (calc_cond)
5492 5334 {
5493 } 5335 char job = '1';
5494 else 5336 Array<Complex> z (2 * nr);
5495 { 5337 Complex *pz = z.fortran_vec ();
5496 if (calc_cond) 5338 Array<double> iz (nr);
5497 { 5339 double *piz = iz.fortran_vec ();
5498 char job = '1'; 5340
5499 Array<Complex> z (2 * nr); 5341 F77_XFCN (zgbcon, ZGBCON,
5500 Complex *pz = z.fortran_vec (); 5342 (F77_CONST_CHAR_ARG2 (&job, 1),
5501 Array<double> iz (nr); 5343 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5502 double *piz = iz.fortran_vec (); 5344 anorm, rcond, pz, piz, err
5503 5345 F77_CHAR_ARG_LEN (1)));
5504 F77_XFCN (zgbcon, ZGBCON, 5346
5505 (F77_CONST_CHAR_ARG2 (&job, 1), 5347 if (err != 0)
5506 nc, n_lower, n_upper, tmp_data, ldm, pipvt, 5348 err = -2;
5507 anorm, rcond, pz, piz, err 5349
5508 F77_CHAR_ARG_LEN (1))); 5350 volatile double rcond_plus_one = rcond + 1.0;
5509 5351
5510 if (f77_exception_encountered) 5352 if (rcond_plus_one == 1.0 || xisnan (rcond))
5511 (*current_liboctave_error_handler) 5353 {
5512 ("unrecoverable error in zgbcon"); 5354 err = -2;
5513 5355
5514 if (err != 0) 5356 if (sing_handler)
5515 err = -2;
5516
5517 volatile double rcond_plus_one = rcond + 1.0;
5518
5519 if (rcond_plus_one == 1.0 || xisnan (rcond))
5520 {
5521 err = -2;
5522
5523 if (sing_handler)
5524 {
5525 sing_handler (rcond);
5526 mattype.mark_as_rectangular ();
5527 }
5528 else
5529 (*current_liboctave_error_handler)
5530 ("matrix singular to machine precision, rcond = %g",
5531 rcond);
5532 }
5533 }
5534 else
5535 rcond = 1.;
5536
5537 if (err == 0)
5538 {
5539 char job = 'N';
5540 volatile octave_idx_type x_nz = b.nnz ();
5541 octave_idx_type b_nc = b.cols ();
5542 retval = SparseComplexMatrix (nr, b_nc, x_nz);
5543 retval.xcidx(0) = 0;
5544 volatile octave_idx_type ii = 0;
5545
5546 OCTAVE_LOCAL_BUFFER (Complex, Bx, nr);
5547
5548 for (volatile octave_idx_type j = 0; j < b_nc; j++)
5549 { 5357 {
5550 for (octave_idx_type i = 0; i < nr; i++) 5358 sing_handler (rcond);
5551 Bx[i] = 0.; 5359 mattype.mark_as_rectangular ();
5552
5553 for (octave_idx_type i = b.cidx(j);
5554 i < b.cidx(j+1); i++)
5555 Bx[b.ridx(i)] = b.data(i);
5556
5557 F77_XFCN (zgbtrs, ZGBTRS,
5558 (F77_CONST_CHAR_ARG2 (&job, 1),
5559 nr, n_lower, n_upper, 1, tmp_data,
5560 ldm, pipvt, Bx, b.rows (), err
5561 F77_CHAR_ARG_LEN (1)));
5562
5563 if (f77_exception_encountered)
5564 {
5565 (*current_liboctave_error_handler)
5566 ("unrecoverable error in dgbtrs");
5567 break;
5568 }
5569
5570 // Count non-zeros in work vector and adjust
5571 // space in retval if needed
5572 octave_idx_type new_nnz = 0;
5573 for (octave_idx_type i = 0; i < nr; i++)
5574 if (Bx[i] != 0.)
5575 new_nnz++;
5576
5577 if (ii + new_nnz > x_nz)
5578 {
5579 // Resize the sparse matrix
5580 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5581 retval.change_capacity (sz);
5582 x_nz = sz;
5583 }
5584
5585 for (octave_idx_type i = 0; i < nr; i++)
5586 if (Bx[i] != 0.)
5587 {
5588 retval.xridx(ii) = i;
5589 retval.xdata(ii++) = Bx[i];
5590 }
5591 retval.xcidx(j+1) = ii;
5592 } 5360 }
5593 5361 else
5594 retval.maybe_compress (); 5362 (*current_liboctave_error_handler)
5595 } 5363 ("matrix singular to machine precision, rcond = %g",
5364 rcond);
5365 }
5366 }
5367 else
5368 rcond = 1.;
5369
5370 if (err == 0)
5371 {
5372 char job = 'N';
5373 volatile octave_idx_type x_nz = b.nnz ();
5374 octave_idx_type b_nc = b.cols ();
5375 retval = SparseComplexMatrix (nr, b_nc, x_nz);
5376 retval.xcidx(0) = 0;
5377 volatile octave_idx_type ii = 0;
5378
5379 OCTAVE_LOCAL_BUFFER (Complex, Bx, nr);
5380
5381 for (volatile octave_idx_type j = 0; j < b_nc; j++)
5382 {
5383 for (octave_idx_type i = 0; i < nr; i++)
5384 Bx[i] = 0.;
5385
5386 for (octave_idx_type i = b.cidx(j);
5387 i < b.cidx(j+1); i++)
5388 Bx[b.ridx(i)] = b.data(i);
5389
5390 F77_XFCN (zgbtrs, ZGBTRS,
5391 (F77_CONST_CHAR_ARG2 (&job, 1),
5392 nr, n_lower, n_upper, 1, tmp_data,
5393 ldm, pipvt, Bx, b.rows (), err
5394 F77_CHAR_ARG_LEN (1)));
5395
5396 // Count non-zeros in work vector and adjust
5397 // space in retval if needed
5398 octave_idx_type new_nnz = 0;
5399 for (octave_idx_type i = 0; i < nr; i++)
5400 if (Bx[i] != 0.)
5401 new_nnz++;
5402
5403 if (ii + new_nnz > x_nz)
5404 {
5405 // Resize the sparse matrix
5406 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5407 retval.change_capacity (sz);
5408 x_nz = sz;
5409 }
5410
5411 for (octave_idx_type i = 0; i < nr; i++)
5412 if (Bx[i] != 0.)
5413 {
5414 retval.xridx(ii) = i;
5415 retval.xdata(ii++) = Bx[i];
5416 }
5417 retval.xcidx(j+1) = ii;
5418 }
5419
5420 retval.maybe_compress ();
5596 } 5421 }
5597 } 5422 }
5598 } 5423 }
5599 else if (typ != MatrixType::Banded_Hermitian) 5424 else if (typ != MatrixType::Banded_Hermitian)
5600 (*current_liboctave_error_handler) ("incorrect matrix type"); 5425 (*current_liboctave_error_handler) ("incorrect matrix type");