This repository has been archived by the owner on Jul 23, 2024. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
03f31a1
commit 4afa62b
Showing
43 changed files
with
9,779 additions
and
7,871 deletions.
There are no files selected for viewing
This file was deleted.
Oops, something went wrong.
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,143 @@ | ||
| /* dpbfa.f -- translated by f2c (version 20190311). | ||
| You must link the resulting object file with libf2c: | ||
| on Microsoft Windows system, link with libf2c.lib; | ||
| on Linux or Unix systems, link with .../path/to/libf2c.a -lm | ||
| or, if you install libf2c.a in a standard place, with -lf2c -lm | ||
| -- in that order, at the end of the command line, as in | ||
| cc *.o -lf2c -lm | ||
| Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., | ||
| http://www.netlib.org/f2c/libf2c.zip | ||
| */ | ||
|
|
||
| #include "f2c.h" | ||
|
|
||
| /* Table of constant values */ | ||
|
|
||
| static integer c__1 = 1; | ||
|
|
||
| /* Subroutine */ int dpbfa_(doublereal *abd, integer *lda, integer *n, | ||
| integer *m, integer *info) | ||
| { | ||
| /* System generated locals */ | ||
| integer abd_dim1, abd_offset, i__1, i__2, i__3; | ||
|
|
||
| /* Builtin functions */ | ||
| double sqrt(doublereal); | ||
|
|
||
| /* Local variables */ | ||
| static integer j, k; | ||
| static doublereal s, t; | ||
| static integer ik, jk, mu; | ||
| extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, | ||
| integer *); | ||
|
|
||
|
|
||
| /* dpbfa factors a double precision symmetric positive definite */ | ||
| /* matrix stored in band form. */ | ||
|
|
||
| /* dpbfa is usually called by dpbco, but it can be called */ | ||
| /* directly with a saving in time if rcond is not needed. */ | ||
|
|
||
| /* on entry */ | ||
|
|
||
| /* abd double precision(lda, n) */ | ||
| /* the matrix to be factored. the columns of the upper */ | ||
| /* triangle are stored in the columns of abd and the */ | ||
| /* diagonals of the upper triangle are stored in the */ | ||
| /* rows of abd . see the comments below for details. */ | ||
|
|
||
| /* lda integer */ | ||
| /* the leading dimension of the array abd . */ | ||
| /* lda must be .ge. m + 1 . */ | ||
|
|
||
| /* n integer */ | ||
| /* the order of the matrix a . */ | ||
|
|
||
| /* m integer */ | ||
| /* the number of diagonals above the main diagonal. */ | ||
| /* 0 .le. m .lt. n . */ | ||
|
|
||
| /* on return */ | ||
|
|
||
| /* abd an upper triangular matrix r , stored in band */ | ||
| /* form, so that a = trans(r)*r . */ | ||
|
|
||
| /* info integer */ | ||
| /* = 0 for normal return. */ | ||
| /* = k if the leading minor of order k is not */ | ||
| /* positive definite. */ | ||
|
|
||
| /* band storage */ | ||
|
|
||
| /* if a is a symmetric positive definite band matrix, */ | ||
| /* the following program segment will set up the input. */ | ||
|
|
||
| /* m = (band width above diagonal) */ | ||
| /* do 20 j = 1, n */ | ||
| /* i1 = max0(1, j-m) */ | ||
| /* do 10 i = i1, j */ | ||
| /* k = i-j+m+1 */ | ||
| /* abd(k,j) = a(i,j) */ | ||
| /* 10 continue */ | ||
| /* 20 continue */ | ||
|
|
||
| /* linpack. this version dated 08/14/78 . */ | ||
| /* cleve moler, university of new mexico, argonne national lab. */ | ||
|
|
||
| /* subroutines and functions */ | ||
|
|
||
| /* blas ddot */ | ||
| /* fortran max0,sqrt */ | ||
|
|
||
| /* internal variables */ | ||
|
|
||
| /* begin block with ...exits to 40 */ | ||
|
|
||
|
|
||
| /* Parameter adjustments */ | ||
| abd_dim1 = *lda; | ||
| abd_offset = 1 + abd_dim1; | ||
| abd -= abd_offset; | ||
|
|
||
| /* Function Body */ | ||
| i__1 = *n; | ||
| for (j = 1; j <= i__1; ++j) { | ||
| *info = j; | ||
| s = 0.; | ||
| ik = *m + 1; | ||
| /* Computing MAX */ | ||
| i__2 = j - *m; | ||
| jk = max(i__2,1); | ||
| /* Computing MAX */ | ||
| i__2 = *m + 2 - j; | ||
| mu = max(i__2,1); | ||
| if (*m < mu) { | ||
| goto L20; | ||
| } | ||
| i__2 = *m; | ||
| for (k = mu; k <= i__2; ++k) { | ||
| i__3 = k - mu; | ||
| t = abd[k + j * abd_dim1] - ddot_(&i__3, &abd[ik + jk * abd_dim1], | ||
| &c__1, &abd[mu + j * abd_dim1], &c__1); | ||
| t /= abd[*m + 1 + jk * abd_dim1]; | ||
| abd[k + j * abd_dim1] = t; | ||
| s += t * t; | ||
| --ik; | ||
| ++jk; | ||
| /* L10: */ | ||
| } | ||
| L20: | ||
| s = abd[*m + 1 + j * abd_dim1] - s; | ||
| /* ......exit */ | ||
| if (s <= 0.) { | ||
| goto L40; | ||
| } | ||
| abd[*m + 1 + j * abd_dim1] = sqrt(s); | ||
| /* L30: */ | ||
| } | ||
| *info = 0; | ||
| L40: | ||
| return 0; | ||
| } /* dpbfa_ */ | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,125 @@ | ||
| /* dpbsl.f -- translated by f2c (version 20190311). | ||
| You must link the resulting object file with libf2c: | ||
| on Microsoft Windows system, link with libf2c.lib; | ||
| on Linux or Unix systems, link with .../path/to/libf2c.a -lm | ||
| or, if you install libf2c.a in a standard place, with -lf2c -lm | ||
| -- in that order, at the end of the command line, as in | ||
| cc *.o -lf2c -lm | ||
| Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., | ||
| http://www.netlib.org/f2c/libf2c.zip | ||
| */ | ||
|
|
||
| #include "f2c.h" | ||
|
|
||
| /* Table of constant values */ | ||
|
|
||
| static integer c__1 = 1; | ||
|
|
||
| /* Subroutine */ int dpbsl_(doublereal *abd, integer *lda, integer *n, | ||
| integer *m, doublereal *b) | ||
| { | ||
| /* System generated locals */ | ||
| integer abd_dim1, abd_offset, i__1, i__2; | ||
|
|
||
| /* Local variables */ | ||
| static integer k; | ||
| static doublereal t; | ||
| static integer kb, la, lb, lm; | ||
| extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, | ||
| integer *); | ||
| extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, | ||
| integer *, doublereal *, integer *); | ||
|
|
||
|
|
||
| /* dpbsl solves the double precision symmetric positive definite */ | ||
| /* band system a*x = b */ | ||
| /* using the factors computed by dpbco or dpbfa. */ | ||
|
|
||
| /* on entry */ | ||
|
|
||
| /* abd double precision(lda, n) */ | ||
| /* the output from dpbco or dpbfa. */ | ||
|
|
||
| /* lda integer */ | ||
| /* the leading dimension of the array abd . */ | ||
|
|
||
| /* n integer */ | ||
| /* the order of the matrix a . */ | ||
|
|
||
| /* m integer */ | ||
| /* the number of diagonals above the main diagonal. */ | ||
|
|
||
| /* b double precision(n) */ | ||
| /* the right hand side vector. */ | ||
|
|
||
| /* on return */ | ||
|
|
||
| /* b the solution vector x . */ | ||
|
|
||
| /* error condition */ | ||
|
|
||
| /* a division by zero will occur if the input factor contains */ | ||
| /* a zero on the diagonal. technically this indicates */ | ||
| /* singularity but it is usually caused by improper subroutine */ | ||
| /* arguments. it will not occur if the subroutines are called */ | ||
| /* correctly and info .eq. 0 . */ | ||
|
|
||
| /* to compute inverse(a) * c where c is a matrix */ | ||
| /* with p columns */ | ||
| /* call dpbco(abd,lda,n,rcond,z,info) */ | ||
| /* if (rcond is too small .or. info .ne. 0) go to ... */ | ||
| /* do 10 j = 1, p */ | ||
| /* call dpbsl(abd,lda,n,c(1,j)) */ | ||
| /* 10 continue */ | ||
|
|
||
| /* linpack. this version dated 08/14/78 . */ | ||
| /* cleve moler, university of new mexico, argonne national lab. */ | ||
|
|
||
| /* subroutines and functions */ | ||
|
|
||
| /* blas daxpy,ddot */ | ||
| /* fortran min0 */ | ||
|
|
||
| /* internal variables */ | ||
|
|
||
|
|
||
| /* solve trans(r)*y = b */ | ||
|
|
||
| /* Parameter adjustments */ | ||
| --b; | ||
| abd_dim1 = *lda; | ||
| abd_offset = 1 + abd_dim1; | ||
| abd -= abd_offset; | ||
|
|
||
| /* Function Body */ | ||
| i__1 = *n; | ||
| for (k = 1; k <= i__1; ++k) { | ||
| /* Computing MIN */ | ||
| i__2 = k - 1; | ||
| lm = min(i__2,*m); | ||
| la = *m + 1 - lm; | ||
| lb = k - lm; | ||
| t = ddot_(&lm, &abd[la + k * abd_dim1], &c__1, &b[lb], &c__1); | ||
| b[k] = (b[k] - t) / abd[*m + 1 + k * abd_dim1]; | ||
| /* L10: */ | ||
| } | ||
|
|
||
| /* solve r*x = y */ | ||
|
|
||
| i__1 = *n; | ||
| for (kb = 1; kb <= i__1; ++kb) { | ||
| k = *n + 1 - kb; | ||
| /* Computing MIN */ | ||
| i__2 = k - 1; | ||
| lm = min(i__2,*m); | ||
| la = *m + 1 - lm; | ||
| lb = k - lm; | ||
| b[k] /= abd[*m + 1 + k * abd_dim1]; | ||
| t = -b[k]; | ||
| daxpy_(&lm, &t, &abd[la + k * abd_dim1], &c__1, &b[lb], &c__1); | ||
| /* L20: */ | ||
| } | ||
| return 0; | ||
| } /* dpbsl_ */ | ||
|
|
Oops, something went wrong.