Actual source code: ilut.c

      1: #define PETSCMAT_DLL

     3:  #include petsc.h
      4: #if !defined(PETSC_AVOID_GNUCOPYRIGHT_CODE)

      6: /* ilut.f -- translated by f2c (version of 25 March 1992  12:58:56).

      8:      The Fortran version of this code was developed by Yousef Saad.
      9:   This code is copyrighted by Yousef Saad with the

     11:                     GNU GENERAL PUBLIC LICENSE
     12:                        Version 2, June 1991

     14:  Copyright (C) 1989, 1991 Free Software Foundation, Inc.
     15:                           675 Mass Ave, Cambridge, MA 02139, USA
     16:  Everyone is permitted to copy and distribute verbatim copies
     17:  of this license document, but changing it is not allowed.

     19:                             Preamble

     21:   The licenses for most software are designed to take away your
     22: freedom to share and change it.  By contrast, the GNU General Public
     23: License is intended to guarantee your freedom to share and change free
     24: software--to make sure the software is free for all its users.  This
     25: General Public License applies to most of the Free Software
     26: Foundation's software and to any other program whose authors commit to
     27: using it.  (Some other Free Software Foundation software is covered by
     28: the GNU Library General Public License instead.)  You can apply it to
     29: your programs, too.

     31:   When we speak of free software, we are referring to freedom, not
     32: price.  Our General Public Licenses are designed to make sure that you
     33: have the freedom to distribute copies of free software (and charge for
     34: this service if you wish), that you receive source code or can get it
     35: if you want it, that you can change the software or use pieces of it
     36: in new free programs; and that you know you can do these things.

     38:   To protect your rights, we need to make restrictions that forbid
     39: anyone to deny you these rights or to ask you to surrender the rights.
     40: These restrictions translate to certain responsibilities for you if you
     41: distribute copies of the software, or if you modify it.

     43:   For example, if you distribute copies of such a program, whether
     44: gratis or for a fee, you must give the recipients all the rights that
     45: you have.  You must make sure that they, too, receive or can get the
     46: source code.  And you must show them these terms so they know their
     47: rights.

     49:   We protect your rights with two steps: (1) copyright the software, and
     50: (2) offer you this license which gives you legal permission to copy,
     51: distribute and/or modify the software.

     53:   Also, for each author's protection and ours, we want to make certain
     54: that everyone understands that there is no warranty for this free
     55: software.  If the software is modified by someone else and passed on, we
     56: want its recipients to know that what they have is not the original, so
     57: that any problems introduced by others will not reflect on the original
     58: authors' reputations.

     60:   Finally, any free program is threatened constantly by software
     61: patents.  We wish to avoid the danger that redistributors of a free
     62: program will individually obtain patent licenses, in effect making the
     63: program proprietary.  To prevent this, we have made it clear that any
     64: patent must be licensed for everyone's free use or not licensed at all.

     66:   The precise terms and conditions for copying, distribution and
     67: modification follow.
     68:
     69:                     GNU GENERAL PUBLIC LICENSE
     70:    TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION

     72:   0. This License applies to any program or other work which contains
     73: a notice placed by the copyright holder saying it may be distributed
     74: under the terms of this General Public License.  The "Program", below,
     75: refers to any such program or work, and a "work based on the Program"
     76: means either the Program or any derivative work under copyright law:
     77: that is to say, a work containing the Program or a portion of it,
     78: either verbatim or with modifications and/or translated into another
     79: language.  (Hereinafter, translation is included without limitation in
     80: the term "modification".)  Each licensee is addressed as "you".

     82: Activities other than copying, distribution and modification are not
     83: covered by this License; they are outside its scope.  The act of
     84: running the Program is not restricted, and the output from the Program
     85: is covered only if its contents constitute a work based on the
     86: Program (independent of having been made by running the Program).
     87: Whether that is true depends on what the Program does.

     89:   1. You may copy and distribute verbatim copies of the Program's
     90: source code as you receive it, in any medium, provided that you
     91: conspicuously and appropriately publish on each copy an appropriate
     92: copyright notice and disclaimer of warranty; keep intact all the
     93: notices that refer to this License and to the absence of any warranty;
     94: and give any other recipients of the Program a copy of this License
     95: along with the Program.

     97: You may charge a fee for the physical act of transferring a copy, and
     98: you may at your option offer warranty protection in exchange for a fee.

    100:   2. You may modify your copy or copies of the Program or any portion
    101: of it, thus forming a work based on the Program, and copy and
    102: distribute such modifications or work under the terms of Section 1
    103: above, provided that you also meet all of these conditions:

    105:     a) You must cause the modified files to carry prominent notices
    106:     stating that you changed the files and the date of any change.

    108:     b) You must cause any work that you distribute or publish, that in
    109:     whole or in part contains or is derived from the Program or any
    110:     part thereof, to be licensed as a whole at no charge to all third
    111:     parties under the terms of this License.

    113:     c) If the modified program normally reads commands interactively
    114:     when run, you must cause it, when started running for such
    115:     interactive use in the most ordinary way, to print or display an
    116:     announcement including an appropriate copyright notice and a
    117:     notice that there is no warranty (or else, saying that you provide
    118:     a warranty) and that users may redistribute the program under
    119:     these conditions, and telling the user how to view a copy of this
    120:     License.  (Exception: if the Program itself is interactive but
    121:     does not normally print such an announcement, your work based on
    122:     the Program is not required to print an announcement.)
    123:
    124: These requirements apply to the modified work as a whole.  If
    125: identifiable sections of that work are not derived from the Program,
    126: and can be reasonably considered independent and separate works in
    127: themselves, then this License, and its terms, do not apply to those
    128: sections when you distribute them as separate works.  But when you
    129: distribute the same sections as part of a whole which is a work based
    130: on the Program, the distribution of the whole must be on the terms of
    131: this License, whose permissions for other licensees extend to the
    132: entire whole, and thus to each and every part regardless of who wrote it.

    134: Thus, it is not the intent of this section to claim rights or contest
    135: your rights to work written entirely by you; rather, the intent is to
    136: exercise the right to control the distribution of derivative or
    137: collective works based on the Program.

    139: In addition, mere aggregation of another work not based on the Program
    140: with the Program (or with a work based on the Program) on a volume of
    141: a storage or distribution medium does not bring the other work under
    142: the scope of this License.

    144:   3. You may copy and distribute the Program (or a work based on it,
    145: under Section 2) in object code or executable form under the terms of
    146: Sections 1 and 2 above provided that you also do one of the following:

    148:     a) Accompany it with the complete corresponding machine-readable
    149:     source code, which must be distributed under the terms of Sections
    150:     1 and 2 above on a medium customarily used for software interchange; or,

    152:     b) Accompany it with a written offer, valid for at least three
    153:     years, to give any third party, for a charge no more than your
    154:     cost of physically performing source distribution, a complete
    155:     machine-readable copy of the corresponding source code, to be
    156:     distributed under the terms of Sections 1 and 2 above on a medium
    157:     customarily used for software interchange; or,

    159:     c) Accompany it with the information you received as to the offer
    160:     to distribute corresponding source code.  (This alternative is
    161:     allowed only for noncommercial distribution and only if you
    162:     received the program in object code or executable form with such
    163:     an offer, in accord with Subsection b above.)

    165: The source code for a work means the preferred form of the work for
    166: making modifications to it.  For an executable work, complete source
    167: code means all the source code for all modules it contains, plus any
    168: associated interface definition files, plus the scripts used to
    169: control compilation and installation of the executable.  However, as a
    170: special exception, the source code distributed need not include
    171: anything that is normally distributed (in either source or binary
    172: form) with the major components (compiler, kernel, and so on) of the
    173: operating system on which the executable runs, unless that component
    174: itself accompanies the executable.

    176: If distribution of executable or object code is made by offering
    177: access to copy from a designated place, then offering equivalent
    178: access to copy the source code from the same place counts as
    179: distribution of the source code, even though third parties are not
    180: compelled to copy the source along with the object code.
    181:
    182:   4. You may not copy, modify, sublicense, or distribute the Program
    183: except as expressly provided under this License.  Any attempt
    184: otherwise to copy, modify, sublicense or distribute the Program is
    185: void, and will automatically terminate your rights under this License.
    186: However, parties who have received copies, or rights, from you under
    187: this License will not have their licenses terminated so long as such
    188: parties remain in full compliance.

    190:   5. You are not required to accept this License, since you have not
    191: signed it.  However, nothing else grants you permission to modify or
    192: distribute the Program or its derivative works.  These actions are
    193: prohibited by law if you do not accept this License.  Therefore, by
    194: modifying or distributing the Program (or any work based on the
    195: Program), you indicate your acceptance of this License to do so, and
    196: all its terms and conditions for copying, distributing or modifying
    197: the Program or works based on it.

    199:   6. Each time you redistribute the Program (or any work based on the
    200: Program), the recipient automatically receives a license from the
    201: original licensor to copy, distribute or modify the Program subject to
    202: these terms and conditions.  You may not impose any further
    203: restrictions on the recipients' exercise of the rights granted herein.
    204: You are not responsible for enforcing compliance by third parties to
    205: this License.

    207:   7. If, as a consequence of a court judgment or allegation of patent
    208: infringement or for any other reason (not limited to patent issues),
    209: conditions are imposed on you (whether by court order, agreement or
    210: otherwise) that contradict the conditions of this License, they do not
    211: excuse you from the conditions of this License.  If you cannot
    212: distribute so as to satisfy simultaneously your obligations under this
    213: License and any other pertinent obligations, then as a consequence you
    214: may not distribute the Program at all.  For example, if a patent
    215: license would not permit royalty-free redistribution of the Program by
    216: all those who receive copies directly or indirectly through you, then
    217: the only way you could satisfy both it and this License would be to
    218: refrain entirely from distribution of the Program.

    220: If any portion of this section is held invalid or unenforceable under
    221: any particular circumstance, the balance of the section is intended to
    222: apply and the section as a whole is intended to apply in other
    223: circumstances.

    225: It is not the purpose of this section to induce you to infringe any
    226: patents or other property right claims or to contest validity of any
    227: such claims; this section has the sole purpose of protecting the
    228: integrity of the free software distribution system, which is
    229: implemented by public license practices.  Many people have made
    230: generous contributions to the wide range of software distributed
    231: through that system in reliance on consistent application of that
    232: system; it is up to the author/donor to decide if he or she is willing
    233: to distribute software through any other system and a licensee cannot
    234: impose that choice.

    236: This section is intended to make thoroughly clear what is believed to
    237: be a consequence of the rest of this License.
    238:
    239:   8. If the distribution and/or use of the Program is restricted in
    240: certain countries either by patents or by copyrighted interfaces, the
    241: original copyright holder who places the Program under this License
    242: may add an explicit geographical distribution limitation excluding
    243: those countries, so that distribution is permitted only in or among
    244: countries not thus excluded.  In such case, this License incorporates
    245: the limitation as if written in the body of this License.

    247:   9. The Free Software Foundation may publish revised and/or new versions
    248: of the General Public License from time to time.  Such new versions will
    249: be similar in spirit to the present version, but may differ in detail to
    250: address new problems or concerns.

    252: Each version is given a distinguishing version number.  If the Program
    253: specifies a version number of this License which applies to it and "any
    254: later version", you have the option of following the terms and conditions
    255: either of that version or of any later version published by the Free
    256: Software Foundation.  If the Program does not specify a version number of
    257: this License, you may choose any version ever published by the Free Software
    258: Foundation.

    260:   10. If you wish to incorporate parts of the Program into other free
    261: programs whose distribution conditions are different, write to the author
    262: to ask for permission.  For software which is copyrighted by the Free
    263: Software Foundation, write to the Free Software Foundation; we sometimes
    264: make exceptions for this.  Our decision will be guided by the two goals
    265: of preserving the free status of all derivatives of our free software and
    266: of promoting the sharing and reuse of software generally.

    268:                             NO WARRANTY

    270:   11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
    271: FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
    272: OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
    273: PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
    274: OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
    275: MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
    276: TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
    277: PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
    278: REPAIR OR CORRECTION.

    280:   12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
    281: WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
    282: REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
    283: INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
    284: OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
    285: TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
    286: YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
    287: PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
    288: POSSIBILITY OF SUCH DAMAGES.

    290:                      END OF TERMS AND CONDITIONS
    291:
    292:         Appendix: How to Apply These Terms to Your New Programs

    294:   If you develop a new program, and you want it to be of the greatest
    295: possible use to the public, the best way to achieve this is to make it
    296: free software which everyone can redistribute and change under these terms.

    298:   To do so, attach the following notices to the program.  It is safest
    299: to attach them to the start of each source file to most effectively
    300: convey the exclusion of warranty; and each file should have at least
    301: the "copyright" line and a pointer to where the full notice is found.

    303:     <one line to give the program's name and a brief idea of what it does.>
    304:     Copyright (C) 19yy  <name of author>

    306:     This program is free software; you can redistribute it and/or modify
    307:     it under the terms of the GNU General Public License as published by
    308:     the Free Software Foundation; either version 2 of the License, or
    309:     (at your option) any later version.

    311:     This program is distributed in the hope that it will be useful,
    312:     but WITHOUT ANY WARRANTY; without even the implied warranty of
    313:     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    314:     GNU General Public License for more details.

    316:     You should have received a copy of the GNU General Public License
    317:     along with this program; if not, write to the Free Software
    318:     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

    320: Also add information on how to contact you by electronic and paper mail.

    322: If the program is interactive, make it output a short notice like this
    323: when it starts in an interactive mode:

    325:     Gnomovision version 69, Copyright (C) 19yy name of author
    326:     Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
    327:     This is free software, and you are welcome to redistribute it
    328:     under certain conditions; type `show c' for details.

    330: The hypothetical commands `show w' and `show c' should show the appropriate
    331: parts of the General Public License.  Of course, the commands you use may
    332: be called something other than `show w' and `show c'; they could even be
    333: mouse-clicks or menu items--whatever suits your program.

    335: You should also get your employer (if you work as a programmer) or your
    336: school, if any, to sign a "copyright disclaimer" for the program, if
    337: necessary.  Here is a sample; alter the names:

    339:   Yoyodyne, Inc., hereby disclaims all copyright interest in the program
    340:   `Gnomovision' (which makes passes at compilers) written by James Hacker.

    342:   <signature of Ty Coon>, 1 April 1989
    343:   Ty Coon, President of Vice

    345: This General Public License does not permit incorporating your program into
    346: proprietary programs.  If your program is a subroutine library, you may
    347: consider it more useful to permit linking proprietary applications with the
    348: library.  If this is what you want to do, use the GNU Library General
    349: Public License instead of this License.

    351: */

    353: static PetscErrorCode SPARSEKIT2qsplit(PetscScalar *a,PetscInt *ind,PetscInt *n,PetscInt *ncut)
    354: {
    355:     /* System generated locals */
    356:     PetscInt i__1;
    357:     PetscScalar d__1;

    359:     /* Local variables */
    360:     PetscInt last,itmp,j,first;
    361:     PetscReal abskey;
    362:     PetscInt mid;
    363:     PetscScalar tmp;

    365: /* -----------------------------------------------------------------------
    366:  */
    367: /*     does a quick-sort split of a real array. */
    368: /*     on input a(1:n). is a real array */
    369: /*     on output a(1:n) is permuted such that its elements satisfy: */

    371: /*     abs(a(i)) .ge. abs(a(ncut)) for i .lt. ncut and */
    372: /*     abs(a(i)) .le. abs(a(ncut)) for i .gt. ncut */

    374: /*    ind(1:n) is an integer array which permuted in the same way as a(*).
    375: */
    376: /* -----------------------------------------------------------------------
    377:  */
    378: /* ----- */
    379:     /* Parameter adjustments */
    380:     --ind;
    381:     --a;

    383:     /* Function Body */
    384:     first = 1;
    385:     last = *n;
    386:     if (*ncut < first || *ncut > last) {
    387:         return 0;
    388:     }

    390: /*     outer loop -- while mid .ne. ncut do */

    392: L1:
    393:     mid = first;
    394:     abskey = (d__1 = a[mid],PetscAbsScalar(d__1));
    395:     i__1 = last;
    396:     for (j = first + 1; j <= i__1; ++j) {
    397:         if ((d__1 = a[j],PetscAbsScalar(d__1)) > abskey) {
    398:             ++mid;
    399: /*     interchange */
    400:             tmp = a[mid];
    401:             itmp = ind[mid];
    402:             a[mid] = a[j];
    403:             ind[mid] = ind[j];
    404:             a[j] = tmp;
    405:             ind[j] = itmp;
    406:         }
    407: /* L2: */
    408:     }

    410: /*     interchange */

    412:     tmp = a[mid];
    413:     a[mid] = a[first];
    414:     a[first] = tmp;

    416:     itmp = ind[mid];
    417:     ind[mid] = ind[first];
    418:     ind[first] = itmp;

    420: /*     test for while loop */

    422:     if (mid == *ncut) {
    423:         return 0;
    424:     }
    425:     if (mid > *ncut) {
    426:         last = mid - 1;
    427:     } else {
    428:         first = mid + 1;
    429:     }
    430:     goto L1;
    431: /* ----------------end-of-qsplit------------------------------------------
    432:  */
    433: /* -----------------------------------------------------------------------
    434:  */
    435: } /* qsplit_ */

    438: /* ---------------------------------------------------------------------- */
    439: PetscErrorCode SPARSEKIT2ilutp(PetscInt *n,PetscScalar *a,PetscInt *ja,PetscInt * ia,PetscInt *lfil,PetscReal droptol,PetscReal *permtol,PetscInt *mbloc,PetscScalar *alu,
    440:         PetscInt *jlu,PetscInt *ju,PetscInt *iwk,PetscScalar *w,PetscInt *jw,  PetscInt *iperm,PetscErrorCode *ierr)
    441: {
    442:     /* System generated locals */
    443:     PetscInt i__1,i__2;
    444:     PetscScalar d__1;

    446:     /* Local variables */
    447:     PetscScalar fact;
    448:     PetscInt lenl,imax,lenu,icut,jpos;
    449:     PetscReal xmax;
    450:     PetscInt jrow;
    451:     PetscReal xmax0;
    452:     PetscInt i,j,k;
    453:     PetscScalar s,t;
    454:     PetscInt j_1,j2;
    455:     PetscReal tnorm,t1;
    456:     PetscInt ii,jj;
    457:     PetscInt ju0,len;
    458:     PetscScalar tmp;

    460: /* -----------------------------------------------------------------------
    461:  */
    462: /*     implicit none */
    463: /* ----------------------------------------------------------------------*
    464:  */
    465: /*       *** ILUTP preconditioner -- ILUT with pivoting  ***            *
    466: */
    467: /*      incomplete LU factorization with dual truncation mechanism      *
    468: */
    469: /* ----------------------------------------------------------------------*
    470:  */
    471: /* author Yousef Saad *Sep 8, 1993 -- Latest revision, August 1996.     *
    472: */
    473: /* ----------------------------------------------------------------------*
    474:  */
    475: /* on entry: */
    476: /* ========== */
    477: /* n       = integer. The dimension of the matrix A. */

    479: /* a,ja,ia = matrix stored in Compressed Sparse Row format. */
    480: /*           ON RETURN THE COLUMNS OF A ARE PERMUTED. SEE BELOW FOR */
    481: /*           DETAILS. */

    483: /* lfil    = integer. The fill-in parameter. Each row of L and each row */

    485: /*           of U will have a maximum of lfil elements (excluding the */
    486: /*           diagonal element). lfil must be .ge. 0. */
    487: /*           ** WARNING: THE MEANING OF LFIL HAS CHANGED WITH RESPECT TO
    488: */
    489: /*           EARLIER VERSIONS. */

    491: /* droptol = real*8. Sets the threshold for dropping small terms in the */

    493: /*           factorization. See below for details on dropping strategy. */

    496: /* lfil    = integer. The fill-in parameter. Each row of L and */
    497: /*           each row of U will have a maximum of lfil elements. */
    498: /*           WARNING: THE MEANING OF LFIL HAS CHANGED WITH RESPECT TO */
    499: /*           EARLIER VERSIONS. */
    500: /*           lfil must be .ge. 0. */

    502: /* permtol = tolerance ratio used to  determne whether or not to permute
    503: */
    504: /*           two columns.  At step i columns i and j are permuted when */

    506: /*                     abs(a(i,j))*permtol .gt. abs(a(i,i)) */

    508: /*           [0 --> never permute; good values 0.1 to 0.01] */

    510: /* mbloc   = if desired, permuting can be done only within the diagonal */

    512: /*           blocks of size mbloc. Useful for PDE problems with several */

    514: /*           degrees of freedom.. If feature not wanted take mbloc=n. */

    517: /* iwk     = integer. The lengths of arrays alu and jlu. If the arrays */
    518: /*           are not big enough to store the ILU factorizations, ilut */
    519: /*           will stop with an error message. */

    521: /* On return: */
    522: /* =========== */

    524: /* alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
    525: */
    526: /*           the L and U factors together. The diagonal (stored in */
    527: /*           alu(1:n)) is inverted. Each i-th row of the alu,jlu matrix
    528: */
    529: /*           contains the i-th row of L (excluding the diagonal entry=1)
    530: */
    531: /*           followed by the i-th row of U. */

    533: /* ju      = integer array of length n containing the pointers to */
    534: /*           the beginning of each row of U in the matrix alu,jlu. */

    536: /* iperm   = contains the permutation arrays. */
    537: /*           iperm(1:n) = old numbers of unknowns */
    538: /*           iperm(n+1:2*n) = reverse permutation = new unknowns. */

    540: /* integer. Error message with the following meaning. */
    541: /*           0    --> successful return. */
    542: /*           ierr .gt. 0  --> zero pivot encountered at step number ierr.
    543: */
    544: /*           -1   --> Error. input matrix may be wrong. */
    545: /*                            (The elimination process has generated a */
    546: /*                            row in L or U whose length is .gt.  n.) */
    547: /*           -2   --> The matrix L overflows the array al. */
    548: /*           -3   --> The matrix U overflows the array alu. */
    549: /*           -4   --> Illegal value for lfil. */
    550: /*           -5   --> zero row encountered. */

    552: /* work arrays: */
    553: /* ============= */
    554: /* jw      = integer work array of length 2*n. */
    555: /* w       = real work array of length n */

    557: /* IMPORTANR NOTE: */
    558: /* -------------- */
    559: /* TO AVOID PERMUTING THE SOLUTION VECTORS ARRAYS FOR EACH LU-SOLVE, */
    560: /* THE MATRIX A IS PERMUTED ON RETURN. [all column indices are */
    561: /* changed]. SIMILARLY FOR THE U MATRIX. */
    562: /* To permute the matrix back to its original state use the loop: */

    564: /*      do k=ia(1), ia(n+1)-1 */
    565: /*         ja(k) = iperm(ja(k)) */
    566: /*      enddo */

    568: /* -----------------------------------------------------------------------
    569:  */
    570: /*     local variables */

    573:     /* Parameter adjustments */
    574:     --iperm;
    575:     --jw;
    576:     --w;
    577:     --ju;
    578:     --jlu;
    579:     --alu;
    580:     --ia;
    581:     --ja;
    582:     --a;

    584:     /* Function Body */
    585:     if (*lfil < 0) {
    586:         goto L998;
    587:     }
    588: /* -----------------------------------------------------------------------
    589:  */
    590: /*     initialize ju0 (points to next element to be added to alu,jlu) */
    591: /*     and pointer array. */
    592: /* -----------------------------------------------------------------------
    593:  */
    594:     ju0 = *n + 2;
    595:     jlu[1] = ju0;

    597: /*  integer PetscReal pointer array. */

    599:     i__1 = *n;
    600:     for (j = 1; j <= i__1; ++j) {
    601:         jw[*n + j] = 0;
    602:         iperm[j] = j;
    603:         iperm[*n + j] = j;
    604: /* L1: */
    605:     }
    606: /* -----------------------------------------------------------------------
    607:  */
    608: /*     beginning of main loop. */
    609: /* -----------------------------------------------------------------------
    610:  */
    611:     i__1 = *n;
    612:     for (ii = 1; ii <= i__1; ++ii) {
    613:         j_1 = ia[ii];
    614:         j2 = ia[ii + 1] - 1;
    615:         tnorm = 0.;
    616:         i__2 = j2;
    617:         for (k = j_1; k <= i__2; ++k) {
    618:             tnorm += (d__1 = a[k], PetscAbsScalar(d__1));
    619: /* L501: */
    620:         }
    621:         if (!tnorm) {
    622:             goto L999;
    623:         }
    624:         tnorm /= j2 - j_1 + 1;

    626: /*     unpack L-part and U-part of row of A in arrays  w  -- */

    628:         lenu = 1;
    629:         lenl = 0;
    630:         jw[ii] = ii;
    631:         w[ii] = (float)0.;
    632:         jw[*n + ii] = ii;

    634:         i__2 = j2;
    635:         for (j = j_1; j <= i__2; ++j) {
    636:             k = iperm[*n + ja[j]];
    637:             t = a[j];
    638:             if (k < ii) {
    639:                 ++lenl;
    640:                 jw[lenl] = k;
    641:                 w[lenl] = t;
    642:                 jw[*n + k] = lenl;
    643:             } else if (k == ii) {
    644:                 w[ii] = t;
    645:             } else {
    646:                 ++lenu;
    647:                 jpos = ii + lenu - 1;
    648:                 jw[jpos] = k;
    649:                 w[jpos] = t;
    650:                 jw[*n + k] = jpos;
    651:             }
    652: /* L170: */
    653:         }
    654:         jj = 0;
    655:         len = 0;

    657: /*     eliminate previous rows */

    659: L150:
    660:         ++jj;
    661:         if (jj > lenl) {
    662:             goto L160;
    663:         }
    664: /* ------------------------------------------------------------------
    665: ----- */
    666: /*     in order to do the elimination in the correct order we must sel
    667: ect */
    668: /*     the smallest column index among jw(k), k=jj+1, ..., lenl. */
    669: /* ------------------------------------------------------------------
    670: ----- */
    671:         jrow = jw[jj];
    672:         k = jj;

    674: /*     determine smallest column index */

    676:         i__2 = lenl;
    677:         for (j = jj + 1; j <= i__2; ++j) {
    678:             if (jw[j] < jrow) {
    679:                 jrow = jw[j];
    680:                 k = j;
    681:             }
    682: /* L151: */
    683:         }

    685:         if (k != jj) {
    686: /*     exchange in jw */
    687:             j = jw[jj];
    688:             jw[jj] = jw[k];
    689:             jw[k] = j;
    690: /*     exchange in jr */
    691:             jw[*n + jrow] = jj;
    692:             jw[*n + j] = k;
    693: /*     exchange in w */
    694:             s = w[jj];
    695:             w[jj] = w[k];
    696:             w[k] = s;
    697:         }

    699: /*     zero out element in row by resetting jw(n+jrow) to zero. */

    701:         jw[*n + jrow] = 0;

    703: /*     get the multiplier for row to be eliminated: jrow */

    705:         fact = w[jj] * alu[jrow];

    707: /*     drop term if small */

    709:         if (PetscAbsScalar(fact) <= droptol) {
    710:             goto L150;
    711:         }

    713: /*     combine current row and row jrow */

    715:         i__2 = jlu[jrow + 1] - 1;
    716:         for (k = ju[jrow]; k <= i__2; ++k) {
    717:             s = fact * alu[k];
    718: /*     new column number */
    719:             j = iperm[*n + jlu[k]];
    720:             jpos = jw[*n + j];
    721:             if (j >= ii) {

    723: /*     dealing with upper part. */

    725:                 if (!jpos) {

    727: /*     this is a fill-in element */

    729:                     ++lenu;
    730:                     i = ii + lenu - 1;
    731:                     if (lenu > *n) {
    732:                         goto L995;
    733:                     }
    734:                     jw[i] = j;
    735:                     jw[*n + j] = i;
    736:                     w[i] = -s;
    737:                 } else {
    738: /*     no fill-in element -- */
    739:                     w[jpos] -= s;
    740:                 }
    741:             } else {

    743: /*     dealing with lower part. */

    745:                 if (!jpos) {

    747: /*     this is a fill-in element */

    749:                     ++lenl;
    750:                     if (lenl > *n) {
    751:                         goto L995;
    752:                     }
    753:                     jw[lenl] = j;
    754:                     jw[*n + j] = lenl;
    755:                     w[lenl] = -s;
    756:                 } else {

    758: /*     this is not a fill-in element */

    760:                     w[jpos] -= s;
    761:                 }
    762:             }
    763: /* L203: */
    764:         }

    766: /*     store this pivot element -- (from left to right -- no danger of
    767:  */
    768: /*     overlap with the working elements in L (pivots). */

    770:         ++len;
    771:         w[len] = fact;
    772:         jw[len] = jrow;
    773:         goto L150;
    774: L160:

    776: /*     reset double-pointer to zero (U-part) */

    778:         i__2 = lenu;
    779:         for (k = 1; k <= i__2; ++k) {
    780:             jw[*n + jw[ii + k - 1]] = 0;
    781: /* L308: */
    782:         }

    784: /*     update L-matrix */

    786:         lenl = len;
    787:         len = PetscMin(lenl,*lfil);

    789: /*     sort by quick-split */

    791:         SPARSEKIT2qsplit(&w[1], &jw[1], &lenl, &len);

    793: /*     store L-part -- in original coordinates .. */

    795:         i__2 = len;
    796:         for (k = 1; k <= i__2; ++k) {
    797:             if (ju0 > *iwk) {
    798:                 goto L996;
    799:             }
    800:             alu[ju0] = w[k];
    801:             jlu[ju0] = iperm[jw[k]];
    802:             ++ju0;
    803: /* L204: */
    804:         }

    806: /*     save pointer to beginning of row ii of U */

    808:         ju[ii] = ju0;

    810: /*     update U-matrix -- first apply dropping strategy */

    812:         len = 0;
    813:         i__2 = lenu - 1;
    814:         for (k = 1; k <= i__2; ++k) {
    815:             if ((d__1 = w[ii + k], PetscAbsScalar(d__1)) > droptol * tnorm) {
    816:                 ++len;
    817:                 w[ii + len] = w[ii + k];
    818:                 jw[ii + len] = jw[ii + k];
    819:             }
    820:         }
    821:         lenu = len + 1;
    822:         len = PetscMin(lenu,*lfil);
    823:         i__2 = lenu - 1;
    824:         SPARSEKIT2qsplit(&w[ii + 1], &jw[ii + 1], &i__2, &len);

    826: /*     determine next pivot -- */

    828:         imax = ii;
    829:         xmax = (d__1 = w[imax], PetscAbsScalar(d__1));
    830:         xmax0 = xmax;
    831:         icut = ii - 1 + *mbloc - (ii - 1) % *mbloc;
    832:         i__2 = ii + len - 1;
    833:         for (k = ii + 1; k <= i__2; ++k) {
    834:             t1 = (d__1 = w[k], PetscAbsScalar(d__1));
    835:             if (t1 > xmax && t1 * *permtol > xmax0 && jw[k] <= icut) {
    836:                 imax = k;
    837:                 xmax = t1;
    838:             }
    839:         }

    841: /*     exchange w's */

    843:         tmp = w[ii];
    844:         w[ii] = w[imax];
    845:         w[imax] = tmp;

    847: /*     update iperm and reverse iperm */

    849:         j = jw[imax];
    850:         i = iperm[ii];
    851:         iperm[ii] = iperm[j];
    852:         iperm[j] = i;

    854: /*     reverse iperm */

    856:         iperm[*n + iperm[ii]] = ii;
    857:         iperm[*n + iperm[j]] = j;
    858: /* ------------------------------------------------------------------
    859: ----- */

    861:         if (len + ju0 > *iwk) {
    862:             goto L997;
    863:         }

    865: /*     copy U-part in original coordinates */

    867:         i__2 = ii + len - 1;
    868:         for (k = ii + 1; k <= i__2; ++k) {
    869:             jlu[ju0] = iperm[jw[k]];
    870:             alu[ju0] = w[k];
    871:             ++ju0;
    872: /* L302: */
    873:         }

    875: /*     store inverse of diagonal element of u */

    877:         if (w[ii] == 0.0) {
    878:             w[ii] = (droptol + 1e-4) * tnorm;
    879:         }
    880:         alu[ii] = 1. / w[ii];

    882: /*     update pointer to beginning of next row of U. */

    884:         jlu[ii + 1] = ju0;
    885: /* ------------------------------------------------------------------
    886: ----- */
    887: /*     end main loop */
    888: /* ------------------------------------------------------------------
    889: ----- */
    890: /* L500: */
    891:     }

    893: /*     permute all column indices of LU ... */

    895:     i__1 = jlu[*n + 1] - 1;
    896:     for (k = jlu[1]; k <= i__1; ++k) {
    897:         jlu[k] = iperm[*n + jlu[k]];
    898:     }

    900: /*     ...and of A */

    902:     i__1 = ia[*n + 1] - 1;
    903:     for (k = ia[1]; k <= i__1; ++k) {
    904:         ja[k] = iperm[*n + ja[k]];
    905:     }

    907:     *0;
    908:     return 0;

    910: /*     incomprehensible error. Matrix must be wrong. */

    912: L995:
    913:     *-1;
    914:     return 0;

    916: /*     insufficient storage in L. */

    918: L996:
    919:     *-2;
    920:     return 0;

    922: /*     insufficient storage in U. */

    924: L997:
    925:     *-3;
    926:     return 0;

    928: /*     illegal lfil entered. */

    930: L998:
    931:     *-4;
    932:     return 0;

    934: /*     zero row encountered */

    936: L999:
    937:     *-5;
    938:     return 0;
    939: /* ----------------end-of-ilutp-------------------------------------------
    940:  */
    941: /* -----------------------------------------------------------------------
    942:  */
    943: } /* ilutp_ */
    944: #endif