From: bangerth Date: Sun, 31 Mar 2013 21:56:39 +0000 (+0000) Subject: Merge from mainline. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9de2bea9b976a5268f51c3d55bb035349f362e30;p=dealii-svn.git Merge from mainline. git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_4@29127 0785d39b-7218-0410-832d-ea1e28bc413d --- 9de2bea9b976a5268f51c3d55bb035349f362e30 diff --cc deal.II/include/deal.II/lac/sparse_direct.h index c4566f101b,195ca99190..3c07c2e8ad --- a/deal.II/include/deal.II/lac/sparse_direct.h +++ b/deal.II/include/deal.II/lac/sparse_direct.h @@@ -1,7 -1,7 +1,7 @@@ //--------------------------------------------------------------------------- // $Id$ // --// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors ++// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@@ -30,998 -30,8 +30,6 @@@ DEAL_II_NAMESPACE_OPEN - /** - * This class provides an interface to the sparse direct solver MA27 - * from the Harwell Subroutine Library. MA27 is a direct solver - * specialized for sparse symmetric indefinite systems of linear - * equations and uses a modified form of Gauss elimination. It is - * included in the Harwell Subroutine - * Library and is written in Fortran. The present class only - * transforms the data stored in SparseMatrix objects into the - * form which is required by the functions resembling MA27, calls - * these Fortran functions, and interprets some of the returned values - * indicating error codes, etc. It also manages allocation of the - * right amount of temporary storage required by these functions. - * - * Note that this class only works if configuration of the deal.II library has - * detected the presence of this solver. Please read the README file on what - * the configure script is looking for and how to provide it. - * - * - *

Interface and Method

- * - * For the meaning of the three functions initialize(), factorize(), - * and solve(), as well as for the method used in MA27, please see the - * documentation - * of these functions. In practice, you will most often call the - * second solve() function, which solves the linear system for a - * given right hand side, but one can as well call the three functions - * separately if, for example, one would like to solve the same matrix - * for several right hand side vectors; the MA27 solver can do this - * efficiently, as it computes a decomposition of the matrix, so that - * subsequent solves only amount to a forward-backward substitution - * which is significantly less costly than the decomposition process. - * - * - *

Parameters to the constructor

- * - * The constructor of this class takes several arguments. The meaning - * is the following: the MA27 functions require the user to allocate - * and pass a certain amount of memory for temporary variables or for - * data to be passed to subsequent functions. The sizes of these - * arrays are denoted by the variables LIW1, LIW2, and LA, - * where LIW1 denotes the size of the IW array in the call to - * MA27A, while LIW2 is the array size in the call to - * MA27B. The documentation of the MA27 functions gives ways to - * obtain estimates for their values, e.g. by evaluating values - * returned by functions called before. However, the documentation - * only states that the values have to be at least as large as - * the estimates, a hint that is not very useful oftentimes (in my - * humble opinion, the lack of dynamic memory allocation mechanism is - * a good reason not to program in Fortran 77 :-). - * - * In our experience, it is often necessary to go beyond the proposed - * values (most often for LA, but also for LIW1). The first - * three parameters of the constructor denote by which factor the - * initial estimates shall be increased. The default values are 1.2 - * (the documentation recommends this value, 1, and 1.5, values which - * have often worked for us. Note that the value of LIW is only - * changed in the second call if the recommended value times - * LIW_factor_2 is larger than the array size already is from the - * call to MA27A; otherwise, LIW_factor_2 is ignored. - * - * If the values thus constructed fail to work, we try to restart the - * called function with larger values until the calls succeed. The - * second triple of values passed to the constructor denotes by which - * factor we shall increase the array sizes. If the increment factors - * are less than or equal to one, then we only try to call the - * respective calls to the functions once and abort by throwing an - * error. Note that the MA27C function writes out an error message - * if the value of LA is too small and gives an indication to - * which size it should be increased. However, most often the - * indicated value is far too small and can not be relied upon. - * - * - *

Note on parallelization

- * - *

Synchronisation

- * - * Due to the use of global variables through COMMON blocks, the calls - * to the sparse direct solver routines are not multithreading-safe, - * i.e. at each time there may only be one call to these functions - * active. You have to synchronise your calls to the functions - * provided by this class using mutexes (see the Threads - * namespace for such classes) to avoid multiple active calls at the - * same time if you use multithreading. Since you may use this class - * in different parts of your program, and may not want to use a - * global variable for locking, this class has a lock as static member - * variable, which may be accessed using the - * get_synchronisation_lock() function. Note however, that this class - * does not perform the synchronisation for you within its member - * functions. The reason is that you will usually want to synchronise - * over the calls to initialize() and factorize(), since there should - * probably not be a call to one of these function with another matrix - * between the calls for one matrix. (The author does not really know - * whether this is true, but it is probably safe to assume that.) - * Since such cross-function synchronisation can only be performed - * from outside, it is left to the user of this class to do so. - * - *

Detached mode

- * - * As an alternative, you can call the function set_detached_mode() - * right after calling the constructor. This lets the program fork, so - * that we now have two programs that communicate via pipes. The - * forked copy of the program then actually replaces itself by a - * program called detached_ma27, that is started in its place - * through the execv system call. Now every time you call one of - * the functions of this class, it relays the data to the other - * program and lets it execute the respective function. The results - * are then transferred back. Since the MA27 functions are only called - * in the detached program, they will now no longer interfere with the - * respective calls to other functions with different data, so no - * synchronisation is necessary any more. - * - * The advantage of this approach is that as many instances of this - * class may be active at any time as you want. This is handy, if your - * programs spens a significant amount of time in them, and you are - * using many threads, for example in a machine with 4 or more - * processors. The disadvantage, of course, is that the data has to - * copied to and from the detached program, which might make things - * slower (though, as we use block writes, this should not be so much - * of a factor). - * - * Since no more synchronisation is necessary, the - * get_synchronisation_lock() returns a reference to a member - * variable when the detached mode is set. Thus, you need not change - * your program: you can still acquire and release the lock as before, - * it will only have no effect now, since different objects of this - * class no longer share the lock, i.e. you will get it always without - * waiting. On the other hand, it will prevent that you call functions - * of this object multiply in parallel at the same time, which is what - * you probably wanted. - * - * - *
Internals of the detached mode
- * - * The program that actually runs the detached solver is called - * detached_ma27, and will show up under this name in the process - * list. It communicates with the main program through a pipe. - * - * Since the solver and the main program are two separated processes, - * the solver program will not be notified if the main program dies, - * for example because it is aborted with Control-C, because an - * exception is raised and not caught, or some other reason. It will - * just not get any new jobs, but will happily wait until the end of - * times. For this reason, the detached solver has a second thread - * running in parallel that simply checks in regular intervals whether - * the main program is still alive, using the ps program. If this - * is no longer the case, the detached solver exits as well. - * - * Since the intervals between two such checks are a couple of second, - * it may happen that the detached solver survives the main program by - * some time. Presently, the check interval is once every 20 - * seconds. After that time, the detached solver should have noticed - * the main programs demise. - * - * @ingroup Solvers Preconditioners - * - * @author Wolfgang Bangerth, 2000, 2001, 2002 - */ - class SparseDirectMA27 : public Subscriptor - { - public: - /** - * Declare type for container size. - */ - typedef types::global_dof_index size_type; - - /** - * Constructor. See the - * documentation of this class - * for the meaning of the - * parameters to this function. - */ - SparseDirectMA27 (const double LIW_factor_1 = 1.2, - const double LIW_factor_2 = 1, - const double LA_factor = 1.5, - const double LIW_increase_factor_1 = 1.2, - const double LIW_increase_factor_2 = 1.2, - const double LA_increase_factor = 1.2, - const bool suppress_output = true); - - /** - * Destructor. - */ - ~SparseDirectMA27 (); - - /** - * Set the detached mode (see the - * general class documentation - * for a description of what this - * is). - * - * This function must not be - * called after initialize() - * (or the two-argument solve() - * function has been called. If - * it is to be called, then only - * right after construction of - * the object, and before first - * use. - */ - void set_detached_mode (); - - /** - * Return whether the detached - * mode is set. - */ - bool detached_mode_set () const; - - /** - * Initialize some data - * structures. This function - * computes symbolically some - * information based on the - * sparsity pattern, but does not - * actually use the values of the - * matrix, so only the sparsity - * pattern has to be passed as - * argument. - */ - void initialize (const SparsityPattern &sparsity_pattern); - - /** - * Actually factorize the - * matrix. This function may be - * called multiple times for - * different matrices, after the - * object of this class has been - * initialized for a certain - * sparsity pattern. You may - * therefore save some computing - * time if you want to invert - * several matrices with the same - * sparsity pattern. However, - * note that the bulk of the - * computing time is actually - * spent in the factorization, so - * this functionality may not - * always be of large benefit. - * - * If the initialization step has - * not been performed yet, then - * the initialize() function is - * called at the beginning of - * this function. - */ - template - void factorize (const SparseMatrix &matrix); - - /** - * Solve for a certain right hand - * side vector. This function may - * be called multiple times for - * different right hand side - * vectors after the matrix has - * been factorized. This yields a - * big saving in computing time, - * since the actual solution is - * fast, compared to the - * factorization of the matrix. - * - * The solution will be returned - * in place of the right hand - * side vector. - * - * If the factorization has not - * happened before, strange - * things will happen. Note that - * we can't actually call the - * factorize() function from - * here if it has not yet been - * called, since we have no - * access to the actual matrix. - */ - template - void solve (Vector &rhs_and_solution) const; - - /** - * Call the three functions - * initialize, factorize and solve - * in that order, i.e. perform - * the whole solution process for - * the given right hand side - * vector. - * - * The solution will be returned - * in place of the right hand - * side vector. - */ - template - void solve (const SparseMatrix &matrix, - Vector &rhs_and_solution); - - /** - * Return an estimate of the - * memory used by this class. - */ - std::size_t memory_consumption () const; - - /** - * Get a reference to the - * synchronisation lock which can - * be used for this class. See - * the general description of - * this class for more - * information. - */ - Threads::Mutex &get_synchronisation_lock () const; - - /** @addtogroup Exceptions - * @{ */ - - /** - * Exception. - */ - DeclException1 (ExcMA27AFailed, - int, - << "The function MA27A failed with an exit code of " << arg1); - /** - * Exception. - */ - DeclException1 (ExcMA27BFailed, - int, - << "The function MA27B failed with an exit code of " << arg1); - /** - * Exception. - */ - DeclException1 (ExcMA27CFailed, - int, - << "The function MA27C failed with an exit code of " << arg1); - /** - * Exception - */ - DeclException0 (ExcInitializeAlreadyCalled); - - /** - * Exception - */ - DeclException0 (ExcFactorizeNotCalled); - - /** - * Exception - */ - DeclException0 (ExcDifferentSparsityPatterns); - - /** - * Exception - */ - DeclException2 (ExcReadError, - int, int, - << "Error while reading in detached mode. Return value " - << "for 'read' was " << arg1 - << ", errno has value " << arg2); - /** - * Exception - */ - DeclException0 (ExcMatrixNotSymmetric); - //@} - private: - /** - * Declare a local type which - * will store the data necessary - * to communicate with a detached - * solver. To avoid adding - * various system include files, - * the actual declaration of this - * class is in the implementation - * file. - */ - struct DetachedModeData; - - /** - * Store in the constructor - * whether the MA27 routines - * shall deliver output to stdout - * or not. - */ - const bool suppress_output; - - /** - * Store whether - * set_detached_mode() has been - * called. - */ - bool detached_mode; - - /** - * Pointer to a structure that - * will hold the data necessary - * to uphold communication with a - * detached solver. - */ - DetachedModeData *detached_mode_data; - - /** - * Store the three values passed - * to the cinstructor. See the - * documentation of this class - * for the meaning of these - * variables. - */ - const double LIW_factor_1; - const double LIW_factor_2; - const double LA_factor; - - /** - * Increase factors in case a - * call to a function fails. - */ - const double LIW_increase_factor_1; - const double LIW_increase_factor_2; - const double LA_increase_factor; - - /** - * Flags storing whether the - * first two functions have - * already been called. - */ - bool initialize_called; - bool factorize_called; - - /** - * Store a pointer to the - * sparsity pattern, to make sure - * that we use the same thing for - * all calls. - */ - SmartPointer sparsity_pattern; - - /** - * Number of nonzero elements in - * the sparsity pattern on and - * above the diagonal. - */ - size_type n_nonzero_elements; - - /** - * Arrays holding row and column - * indices. - */ - std::vector row_numbers; - std::vector column_numbers; - - /** - * Array to hold the matrix - * elements, and later the - * elements of the factors. - */ - std::vector A; - - /** - * Length of the A array. - */ - size_type LA; - - /** - * Scratch arrays and variables - * used by the MA27 functions. We - * keep to the names introduced - * in the documentation of these - * functions, in all uppercase - * letters as is usual in - * Fortran. - */ - size_type LIW; - std::vector IW; - std::vector IKEEP; - std::vector IW1; - - size_type NSTEPS; - size_type MAXFRT; - - /** - * Two values that live inside a - * COMMON block of the Fortran - * code and are mirrored at these - * locations. They are used to - * transport information about - * the required length of arrays - * from the Fortran functions to - * the outside world. - */ - size_type NRLNEC; - size_type NIRNEC; - - /** - * Flag indicating the level of - * output desired and returning - * error values if error occured. - */ - int IFLAG; - - /** - * Mutexes for synchronising access - * to this class. - */ - static Threads::Mutex static_synchronisation_lock; - mutable Threads::Mutex non_static_synchronisation_lock; - - /** - * Fill the A array from the - * symmetric part of the given - * matrix. - */ - template - void fill_A (const SparseMatrix &matrix); - - /** - * Call the respective function - * with the given args, either - * locally or remote. - */ - void call_ma27ad (const size_type *N, - const size_type *NZ, - const size_type *IRN, - const size_type *ICN, - size_type *IW, - const size_type *LIW, - size_type *IKEEP, - size_type *IW1, - size_type *NSTEPS, - int *IFLAG); - - /** - * Call the respective function - * with the given args, either - * locally or remote. - */ - void call_ma27bd (const size_type *N, - const size_type *NZ, - const size_type *IRN, - const size_type *ICN, - double *A, - const size_type *LA, - size_type *IW, - const size_type *LIW, - const size_type *IKEEP, - const size_type *NSTEPS, - size_type *MAXFRT, - size_type *IW1, - int *IFLAG); - - /** - * Call the respective function - * with the given args, either - * locally or remote. - */ - void call_ma27cd (const size_type *N, - const double *A, - const size_type *LA, - const size_type *IW, - const size_type *LIW, - const size_type *MAXFRT, - double *RHS, - const size_type *IW1, - const size_type *NSTEPS) const; - - /** - * Call the respective function - * with the given args, either - * locally or remote. - */ - void call_ma27x1 (size_type *NRLNEC); - - /** - * Call the respective function - * with the given args, either - * locally or remote. - */ - void call_ma27x2 (size_type *NIRNEC); - - /** - * Call the respective function - * with the given args, either - * locally or remote. - */ - void call_ma27x3 (const size_type *LP); - }; - - - - - - - /** - * This class provides an interface to the sparse direct solver MA47 - * from the Harwell Subroutine Library. MA47 is a direct solver - * specialized for sparse symmetric indefinite systems of linear - * equations and uses a frontal elimination method. It is included in - * the Harwell - * Subroutine Library and is written in Fortran. The present class - * only transforms the data stored in SparseMatrix objects into - * the form which is required by the functions resembling MA47, calls - * these Fortran functions, and interprets some of the returned values - * indicating error codes, etc. It also manages allocation of the - * right amount of temporary storage required by these functions. - * - * Note that this class only works if configuration of the deal.II library has - * detected the presence of this solver. Please read the README file on what - * the configure script is looking for and how to provide it. - * - * - *

Interface and Method

- * - * For the meaning of the three functions initialize(), factorize(), - * and solve(), as well as for the method used in MA47, please see the - * documentation - * of these functions. In practice, one will most often call the - * second solve() function, which solves the linear system for a given - * right hand side, but one can as well call the three functions - * separately if, for example, one would like to solve the same matrix - * for several right hand side vectors; the MA47 solver can do this - * efficiently, as it computes a decomposition of the matrix, so that - * subsequent solves only amount to a forward-backward substitution - * which is significantly less costly than the decomposition process. - * - * - *

Parameters to the constructor

- * - * The constructor of this class takes several arguments. Their - * meaning is equivalent to those of the constructor of the - * SparseDirectMA27 class; see there for more information. - * - * - *

Note on parallelization

- * - * Due to the use of global variables through COMMON blocks, the calls - * to the sparse direct solver routines is not multithreading-capable, - * i.e. at each time there may only be one call to these functions - * active. You have to synchronise your calls to the functions - * provided by this class using mutexes (see the Threads - * namespace for such classes) to avoid multiple active calls at the - * same time if you use multithreading. Since you may use this class - * in different parts of your program, and may not want to use a - * global variable for locking, this class has a lock as static member - * variable, which may be accessed using the - * get_synchronisation_lock() function. Note however, that this class - * does not perform the synchronisation for you within its member - * functions. The reason is that you will usually want to synchronise - * over the calls to initialize() and factorize(), since there should - * probably not be a call to one of these function with another matrix - * between the calls for one matrix. (The author does not really know - * whether this is true, but it is probably safe to assume that.) - * Since such cross-function synchronisation can only be performed - * from outside, it is left to the user of this class to do so. - * - * A detached mode as for MA27 has not yet been implemented for this - * class. - * - * - * @ingroup Solvers Preconditioners - * - * @author Wolfgang Bangerth, 2000, 2001 - */ - class SparseDirectMA47 : public Subscriptor - { - public: - /** - * Declare type for container size. - */ - typedef types::global_dof_index size_type; - - /** - * Constructor. See the - * documentation of this class - * for the meaning of the - * parameters to this function. - * - * This function already calls - * the initialization function - * MA47ID to set up some - * values. - */ - SparseDirectMA47 (const double LIW_factor_1 = 1.4, - const double LIW_factor_2 = 1, - const double LA_factor = 3, - const double LIW_increase_factor_1 = 1.2, - const double LIW_increase_factor_2 = 1.2, - const double LA_increase_factor = 1.2, - const bool suppress_output = true); - - /** - * Initialize some data - * structures. This function - * computes symbolically some - * information based on the - * sparsity pattern, but does not - * actually use the values of the - * matrix, so only the sparsity - * pattern has to be passed as - * argument. - * - * Since the MA47 solver requires - * us to omit zero-entries in the - * matrix (even if they are in - * the sparsity pattern), we have - * to actually use the matrix - * here, as opposed to the MA27 - * solver that only required the - * sparsity pattern. - */ - void initialize (const SparseMatrix &matrix); - - /** - * Actually factorize the - * matrix. Unlike for the MA27 - * solver, this function may not - * be called multiple times for - * different matrices, since we - * have eliminated entries from - * the sparsity pattern where - * matrix entries happen to be - * zero. Since this is likely to - * change between matrices - * although they have the same - * sparsity pattern. - * - * If the initialization step has - * not been performed yet, then - * the initialize() function is - * called at the beginning of - * this function. - */ - void factorize (const SparseMatrix &matrix); - - /** - * Solve for a certain right hand - * side vector. This function may - * be called multiple times for - * different right hand side - * vectors after the matrix has - * been factorized. This yields a - * big saving in computing time, - * since the actual solution is - * fast, compared to the - * factorization of the matrix. - * - * The solution will be returned - * in place of the right hand - * side vector. - * - * If the factorization has not - * happened before, strange - * things will happen. Note that - * we can't actually call the - * factorize() function from - * here if it has not yet been - * called, since we have no - * access to the actual matrix. - */ - void solve (Vector &rhs_and_solution); - - /** - * Call the three functions - * initialize, factorize and - * solve - * in that order, i.e. perform - * the whole solution process for - * the given right hand side - * vector. - * - * The solution will be returned - * in place of the right hand - * side vector. - */ - void solve (const SparseMatrix &matrix, - Vector &rhs_and_solution); - - /** - * Return an estimate of the - * memory used by this class. - */ - std::size_t memory_consumption () const; - - /** - * Get a reference to the - * synchronisation lock which can - * be used for this class. See - * the general description of - * this class for more - * information. - */ - Threads::Mutex &get_synchronisation_lock () const; - - /** @addtogroup Exceptions - * @{ */ - - /** - * Exception. - */ - DeclException1 (ExcMA47AFailed, - int, - << "The function MA47A failed with an exit code of " << arg1); - /** - * Exception. - */ - DeclException1 (ExcMA47BFailed, - int, - << "The function MA47B failed with an exit code of " << arg1); - /** - * Exception. - */ - DeclException1 (ExcMA47CFailed, - int, - << "The function MA47C failed with an exit code of " << arg1); - /** - * Exception - */ - DeclException0 (ExcInitializeAlreadyCalled); - - /** - * Exception - */ - DeclException0 (ExcFactorizeNotCalled); - - /** - * Exception - */ - DeclException0 (ExcCantFactorizeAgain); - - /** - * Exception - */ - DeclException0 (ExcDifferentMatrices); - /** - * Exception - */ - DeclException0 (ExcMatrixNotSymmetric); - //@} - private: - /** - * Store in the constructor - * whether the MA47 routines - * shall deliver output to stdout - * or not. - */ - const bool suppress_output; - - /** - * Store the three values passed - * to the cinstructor. See the - * documentation of this class - * for the meaning of these - * variables. - */ - const double LIW_factor_1; - const double LIW_factor_2; - const double LA_factor; - - /** - * Increase factors in case a - * call to a function fails. - */ - const double LIW_increase_factor_1; - const double LIW_increase_factor_2; - const double LA_increase_factor; - - /** - * Flags storing whether the - * first two functions have - * already been called. - */ - bool initialize_called; - bool factorize_called; - - /** - * Store a pointer to the matrix, - * to make sure that we use the - * same thing for all calls. - */ - SmartPointer,SparseDirectMA47> matrix; - - /** - * Number of nonzero elements in - * the sparsity pattern on and - * above the diagonal. - */ - size_type n_nonzero_elements; - - /** - * Control values set by MA47ID. - */ - double CNTL[2]; - unsigned int ICNTL[7]; - - /** - * Info field filled by the MA47 - * functions and (partially) used - * for subsequent MA47 calls. - */ - int INFO[24]; - - /** - * Arrays holding row and column - * indices. - */ - std::vector row_numbers; - std::vector column_numbers; - - /** - * Array to hold the matrix - * elements, and later the - * elements of the factors. - */ - std::vector A; - - /** - * Length of the A array. - */ - size_type LA; - - /** - * Scratch arrays and variables - * used by the MA47 functions. We - * keep to the names introduced - * in the documentation of these - * functions, in all uppercase - * letters as is usual in - * Fortran. - */ - size_type LIW; - std::vector IW; - std::vector KEEP; - std::vector IW1; - - /** - * Mutex for synchronising access - * to this class. - */ - static Threads::Mutex synchronisation_lock; - - /** - * Fill the A array from the - * symmetric part of the given - * matrix. - */ - void fill_A (const SparseMatrix &matrix); - - /** - * Call the ma47id function - * with the given args. - */ - void call_ma47id (double *CNTL, - unsigned int *ICNTL); - - /** - * Call the ma47ad function - * with the given args. - */ - void call_ma47ad (const size_type *n_rows, - const size_type *n_nonzero_elements, - size_type *row_numbers, - size_type *column_numbers, - size_type *IW, - const size_type *LIW, - size_type *KEEP, - const unsigned int *ICNTL, - int *INFO); - - /** - * Call the ma47bd function - * with the given args. - */ - void call_ma47bd (const size_type *n_rows, - const size_type *n_nonzero_elements, - const size_type *column_numbers, - double *A, - const size_type *LA, - size_type *IW, - const size_type *LIW, - const size_type *KEEP, - const double *CNTL, - const unsigned int *ICNTL, - size_type *IW1, - int *INFO); - - /** - * Call the ma47bd function - * with the given args. - */ - void call_ma47cd (const size_type *n_rows, - const double *A, - const size_type *LA, - const size_type *IW, - const size_type *LIW, - double *rhs_and_solution, - size_type *IW1, - const unsigned int *ICNTL); - }; - - -- -- /** * This class provides an interface to the sparse direct solver * UMFPACK (see matrix. + * This function initializes a MUMPS instance and hands over the system's + * matrix matrix. */ template void initialize_matrix (const Matrix &matrix); diff --cc deal.II/source/lac/sparse_direct.cc index 1e8bd2dbff,1328eb20ce..f7625897db --- a/deal.II/source/lac/sparse_direct.cc +++ b/deal.II/source/lac/sparse_direct.cc @@@ -1550,38 -101,26 +101,26 @@@ voi SparseDirectUMFPACK:: sort_arrays (const SparseMatrix &matrix) { - // do the copying around of entries - // so that the diagonal entry is in the - // right place. note that this is easy to - // detect: since all entries apart from the - // diagonal entry are sorted, we know that - // the diagonal entry is in the wrong place - // if and only if its column index is - // larger than the column index of the - // second entry in a row + // do the copying around of entries so that the diagonal entry is in the + // right place. note that this is easy to detect: since all entries apart + // from the diagonal entry are sorted, we know that the diagonal entry is + // in the wrong place if and only if its column index is larger than the + // column index of the second entry in a row // // ignore rows with only one or no entry - for (unsigned int row=0; row Ai[cursor+1])) @@@ -1621,25 -160,17 +160,17 @@@ voi SparseDirectUMFPACK:: sort_arrays (const BlockSparseMatrix &matrix) { - // the case for block matrices is a - // bit more difficult, since all we - // know is that *within each - // block*, the diagonal of that - // block may come first. however, - // that means that there may be as - // many entries per row in the - // wrong place as there are block - // columns. we can do the same - // thing as above, but we have to - // do it multiple times + // the case for block matrices is a bit more difficult, since all we know + // is that *within each block*, the diagonal of that block may come + // first. however, that means that there may be as many entries per row + // in the wrong place as there are block columns. we can do the same + // thing as above, but we have to do it multiple times - for (unsigned int row=0; row(Ap.back()) == Ai.size(), + Assert (static_cast(Ap.back()) == Ai.size(), ExcInternalError()); - // then copy over matrix - // elements. note that for sparse - // matrices, iterators are sorted - // so that they traverse each row - // from start to end before moving - // on to the next row. however, - // this isn't true for block - // matrices, so we have to do a bit - // of book keeping + // then copy over matrix elements. note that for sparse matrices, + // iterators are sorted so that they traverse each row from start to end + // before moving on to the next row. however, this isn't true for block + // matrices, so we have to do a bit of book keeping { - // have an array that for each - // row points to the first entry - // not yet written to + // have an array that for each row points to the first entry not yet + // written to std::vector row_pointers = Ap; - // loop over the elements of the matrix row by row, as suggested - // in the documentation of the sparse matrix iterator class + // loop over the elements of the matrix row by row, as suggested in the + // documentation of the sparse matrix iterator class - for (unsigned int row = 0; row < matrix.m(); ++row) + for (size_type row = 0; row < matrix.m(); ++row) { for (typename Matrix::const_iterator p=matrix.begin(row); p!=matrix.end(row); ++p) @@@ -1749,9 -260,8 +260,8 @@@ } } - // at the end, we should have - // written all rows completely + // at the end, we should have written all rows completely - for (unsigned int i=0; i