Alexandre Bardakoff

Bruno Bachelet

Timothy Blattner

Walid Keyrouz

Gerson C Kroiz

Loïc Yon

Hedgehog: Understandable Scheduler-Free Heterogeneous Asynchronous Multithreaded Data-Flow Graphs

Getting performance on high-end heterogeneous nodes is challenging. This is due to the large semantic gap between a computation's specification-possibly mathematical formulas or an abstract sequential algorithm-and its parallel implementation; this gap obscures the program's parallel structures and how it gains or loses performance. We present Hedgehog, a library aimed at coarse-grain parallelism. It explicitly embeds a data-flow graph in a program and uses this graph at runtime to drive the program's execution so it takes advantage of hardware parallelism (multicore CPUs and multiple accelerators). Hedgehog has asynchronicity built in. It statically binds individual threads to graph nodes, which are ready to fire when any of their inputs are available. This allows Hedgehog to avoid using a global scheduler and the loss of performance associated with global synchronizations and managing of thread pools. Hedgehog provides a separation of concerns and distinguishes between compute and state maintenance tasks. Its API reflects this separation and allows a developer to gain a better understanding of performance when executing the graph. Hedgehog is implemented as a C ++ 17 headers-only library. One feature of the framework is its low overhead; it transfers control of data between two nodes in ≈ 1 µs. This low overhead combines with Hedgehog's API to provide essentially cost-free profiling of the graph, thereby enabling experimentation for performance, which enhances a developer's insight into a program's performance.

Hedgehog's asynchronous data-flow graph supports a data streaming programming model both within and between graphs. We demonstrate the effectiveness of this approach by highlighting the performance of streaming implementations of two numerical linear algebra routines, which are comparable to existing libraries: matrix multiplication achieves >95 % of the theoretical peak of 4 GPUs; LU decomposition with partial pivoting starts streaming partial final result blocks 40× earlier than waiting for the full result. The relative ease and understandability of obtaining performance with Hedgehog promises to enable non-specialists to target performance on high-end single nodes.

Introduction

Parallel programs are increasing in complexity as they target massively parallel platforms (e.g., 2 × 64 or more CPU cores and multiple GPUs). Obtaining performance on such platforms is challenging because of the large semantic gap between a computation's specification and its implementation. The specification is often as simple as a set of mathematical formulas or an abstract sequential algorithm. By contrast, a parallel implementation must address several coupled issues beyond what a sequential implementation addresses: concurrent computations, multiple memory resources, race conditions and deadlocks, and data motion costs. One approach to bridge this gap is to use a framework that simplifies application development and, equally important, exposes abstractions that address parallelism and performance as first-class concerns. Furthermore, these abstractions should represent parallel constructs and make it easier to instrument and reason about an application's performance thereby allowing developers to gain deeper insight. This paper presents Hedgehog, a general-purpose performance-oriented library aimed at coarsegrain parallelism on single high-end heterogeneous compute nodes. A Hedgehog program contains an explicit representation of a static data-flow graph. The graph executes in a pure asynchronous data-driven mode without a global scheduler and statically binds threads to persistent tasks in the graph. Hedgehog provides a separation of concerns and distinguishes between compute tasks (i.e., compute-bound kernels) and state manager tasks. State managers represent and manage local states between the compute tasks that they connect to. Hedgehog also provides a memory manager tool to control the use of memory resources within a task.

Hedgehog transfers control of data efficiently between two tasks (≈ 1 µs). Its explicit data-flow representation and the graph's cost-free profiling allow the developer to quickly identify how the overall computation is carried out. In addition, it encourages a developer to experiment with the data-flow graph itself and with how to customize the degree of parallelism of compute tasks.

The rest of this paper is organized as follows. Section 2 presents solution requirements. Section 3 reviews existing approaches and frameworks. Section 4 discusses Hedgehog and provides a high-level view of its implementation. Section 5 examines two linear algebra examples using Hedgehog and highlights obtained performances. Section 6 concludes and Section 7 outlines future work.

Solution Approach

Developing an application to scale on a high-end heterogeneous node is challenging. These nodes have a high degree of hardware parallelism with high core-count CPUs (currently 2× 64-core CPUs, 128-core CPUs soon), multiple accelerators and GPUs (up to 4 or 8 per node), and multi-instance GPUs (7× per Nvidia A100 GPU). Furthermore, these nodes have become commodity and, as such, are now accessible to a broad community that is much disconnected from the traditional HPC community and lacks its deep knowledge, which is needed to take advantage of these nodes. An enabling solution for this nascent HPC community should compose optimized compute kernels and efficient data transfers between memory domains (CPU and GPU), express the locality of data, overlap data motion and computations, and fully utilize available hardware. Equally important, this solution should have an explicit and understandable program and execution model that is accessible to both expert and non-expert developers so they can improve their code's performance. This model will allow developers to operate at higher levels of abstractions and will also make it easier for them to adapt their software designs to future architectures as hardware evolves.

Existing approaches can exploit the compute power available in high-end nodes, but do so with an implicit execution model that requires the use of external profiling tools to reveal the impact of software design decisions on performance. This makes it challenging for developers to reason about optimization strategies. Furthermore, these external tools may have substantial overhead; this may make it even more challenging for developers to pinpoint performance bottlenecks.

Our solution, Hedgehog, is based on a static and explicit data-flow graph that operates without a global scheduler and is aimed at coarse-grain parallelism. The execution model obtains performance via data pipelining strategies, which works best with streaming data. Using this approach, the developer can effectively overlap data motion and computations to keep the hardware busy. The model is accompanied by tools that assist with memory management to express data locality in multiple memory domains (CPU and GPUs), thereby addressing memory constraints on a heterogeneous node. We use the explicit representation throughout the execution and profiling of the program; this leads to cost-free visual feedback. The representation and feedback make it easier to reason about the graph's execution and encourage the developer to experiment with the graph so as to optimize its performance. This experimentation can lead to quickly identifying the performance critical path in an implementation's graph. We designate this approach as portable designs for performance because the underlying data-flow graph can be easily tuned to different high-end nodes. These decisions can be analyzed in conjunction with the graph profiling to identify optimal configurations. This iterative optimization approach also applies to improving tasks to take better advantage of hardware, also called experimentation for performance.

Background/State of the Art

Multiple models exist for developing parallel applications that maximize the utilization of available hardware and can target future heterogeneous nodes. The three most commonly used approaches are parallel libraries, language extensions, and task-based libraries.

Parallelized libraries, such as OpenBLAS [START_REF] Wang | AUGEM: Automatically generate high performance Dense Linear Algebra kernels on x86 CPUs[END_REF], OpenCV [START_REF] Culjak | A brief introduction to OpenCV[END_REF], or FFTW [START_REF] Frigo | The Design and Implementation of FFTW3[END_REF], implement parallelism within each function call and effectively establish a synchronization barrier at the end of each function invocation in the parallelized library. Hedgehog encourages the use of parallelized libraries in single-threaded mode within its tasks for their optimized implementations that use hardware vector instructions for example.

Directive-based approaches, such as OpenMP [START_REF]OpenMP API Specification: Version 5[END_REF], OpenACC [START_REF] Herdman | Achieving Portability and Performance through OpenACC[END_REF], and OmpSs-2 [START_REF]The OmpSs-2 Programming Model[END_REF] use mostly pragmas or codelets. These approaches focus on using loop parallelism to obtain performance. They offer many options to customize the pragmas, which can become very complex when handling loops with multiple dependencies and synchronization points. Therefore, to be used correctly, they often require power users with a deep understanding of the hardware and have knowledge about such parallel programming techniques. The data-flow model is capable of fully utilizing a high-end node by orchestrating coarse-grain task parallelism. It also allows an end-user to invoke parallel directive-based kernels from existing libraries in a task.

Task-based approaches exist under numerous forms with different properties. The majority of task-based libraries operate with two prevalent traits: (1) the graph's representation as a DAG (Directed Acyclic Graph) and ( 2) the usage of a pool of threads. The programmer describes the DAG by specifying dependencies between tasks, which task is applicable to either a processor, co-processor, or both. Efficiently binding these tasks to threads in the thread pool, which also selects either the processor or co-processor as the compute resource, is non-trivial. There is no perfect dynamic matching algorithm, which is considered a NP-complete problem, but optimizers can be used to improve the matching. Furthermore, to speed up the overall computation, many task-based approaches add work-stealing or work-balancing techniques at the cost of increasing overhead. These steps define how the library intrinsically works, and are hidden from the algorithm developer in the interest of simplifying parallel programming.

Existing task-based libraries, such as HPX [START_REF] Kaiser | STEllAR-GROUP/hpx: HPX V1.4.1: The C++ Standards Library for Parallelism and Concurrency[END_REF], Legion [START_REF] Bauer | Legion: Expressing Locality and Independence with Logical Regions[END_REF], StarPU [START_REF] Augonnet | StarPU: A Unified Platform for Task Scheduling on Heterogeneous Multicore Architectures[END_REF] or Charm++ [START_REF] Laxmitant | The CHARM Parallel Programming Language and System: Part I -Description of Language Features[END_REF], target performance at any scale from fine-grained parallelism to distributed parallel computing. HPX uses an API compatible with the C ++ 14 standard and is based on message-driven computations and asynchronous calls. Legion bases its computation representation on the decomposition of its data. The logical regions express locality and dependence of program data and tasks that act on the logical regions. StarPU uses a codelet approach. Charm++ is a message-driven library that depends on asynchronous calls. Hedgehog focuses on maximizing hardware utilization on a single node, while maintaining low latency, from 1 µs to 10 µs for a large variety of tasks. Tasks are represented as C ++ classes that embed useful information for the computation (class attributes, static variables, usage of input/output streams) as opposed to using codelets.

Specialized libraries, such as Uintah [START_REF] Meng | Uintah Hybrid Task-Based Parallelism Algorithm[END_REF], Kokkos [START_REF] Edwards | Kokkos[END_REF], or Halide [START_REF] Ragan-Kelley | Halide: A Language and Compiler for Optimizing Parallelism, Locality, and Recomputation in Image Processing Pipelines[END_REF], have been developed with key design decisions for specific application classes. Uintah is based on adaptive mesh refinement with a runtime system aimed at solving partial differential equations in large scale simulations. Kokkos, a library used for manycore parallelism on clusters with MPI, targets HPC applications. It also comes with a co-library, kokkos-kernel, specialized for linear algebra. Halide, a language embedded in C ++ , targets image processing algorithms. Hedgehog has no specialization, but benefits most from an effective coarse-grained data decomposition that feeds sufficient data into the data-flow graph.

Intel Threading Building Blocks [START_REF] Kukanov | The Foundations for Scalable Multicore Software in Intel Threading Building Blocks[END_REF] is a template library for parallel programming. It proposes an algorithm as skeletons along with parallel containers to achieve performance on general algorithms. An algorithm is automatically represented internally as a graph with the library in charge of managing the algorithm's graph from creation to tear-down with work stealing for load-balancing on the node. Hedgehog uses an explicit graph to construct the algorithm and maintains it during execution. No scheduling or balancing are used because Hedgehog relies only on its data-flow representation.

Hedgehog evolved from HTGS [START_REF] Blattner | A hybrid task graph scheduler for high performance image processing workflows[END_REF]; they are a part of task-based approach with some major differences explained in Section 4. Hedgehog is at least as efficient as HTGS, according to a study conducted in Section 5.1. There is no loss of features; instead some features have been added such as multiple inputs, broadcast, and an internal re-architecture transforming the HTGS Bookkeeper into the State Manager.

Hedgehog

Hedgehog is a software infrastructure that presents an algorithm at a high level of abstraction and helps developers reason about their applications. The library is header-only and implemented using the C ++ 17 standard, and requires no other external dependencies.

Algorithms are formulated into a data-flow graph using the Hedgehog library. Nodes in the graph are persistent entities that accept and produce data. Edges connect nodes using queues that store data.

Nodes in Hedgehog have been designed for non-overlapping usage to achieve a separation of concerns. A task is a node that carries out computations, where each task is statically bound to a thread. The state manager is a specialized task that embeds a shareable thread-safe state object to locally manage and do synchronization on the computation flow. This separation of concerns is a cornerstone of Hedgehog's design and facilitates the programmability of the library. The nodes and edges operate without any added scheduler and, alongside the threads attached to each task, formulate an asynchronous data pipeline, which allows Hedgehog to overlap computation, data motion, and I/O. It leverages streaming execution using data decomposition to maximize the system utilization.

The graph is also a node and can be embedded in another graph; this structure enables composition and code sharing. The Hedgehog graph is static so there is no back-end reorganization or expansion of the graph.

To present Hedgehog and its model, we first present the data-flow model in Section 4.1, with the different kinds of nodes in Section 4.2 and their interfaces in Section 4.3. Section 4.4 uncovers our threading model and how the computation is conducted without a scheduler. Section 4.5 and Section 4.6 reveal tools embedded in Hedgehog such as the memory manager, and our costfree visual feedback. Section 4.7 exposes how Hedgehog enforces type compatibility with C ++ metaprogramming techniques.

To illustrate the various components of the Hedgehog library, we present a CPU-based implementation of matrix multiplication. In Section 5, we present the results of a variant of this implementation that targets multi-GPU nodes to showcase the performance and low overhead of Hedgehog's operation. To be clear, Hedgehog is a general purpose library; its programming model has been used in a variety of applications such as image processing and signal processing.

In this example of matrix multiplication, we will create a BLAS-like routine, similar to GEneral Matrix Multiplication (GEMM ), C = A × B + C with matrices A n×m , B m×p , and C n×p . These matrices are decomposed into blocks: C is decomposed into N × P blocks along the rows and columns of the matrix, respectively. In order to access a block at row r and column c we use the notation C r,c . The computation is achieved by iterating for each block in

C, C r,c = M i=1 (A r,i × B i,c
), where M is the number of blocks along the shared dimension. This example as well as others are available on github [START_REF] Bardakoff | Hedgehog Tutorials[END_REF]. These tutorials demonstrate the API.

Structural model: data-flow graph

The first step to programming with Hedgehog is to understand how data flows within an algorithm. The best way to approach this step is to formulate the algorithm as a data-flow graph, which is a representation for computation using a directed graph. The graph's nodes are the computation actors and the graph's edges are the directed information flow. The graph has one entry point and one exit point.

Every task executes based on the presence of data in its input queues, which results in a pipeline that executes concurrently. This threading model is described with more detail in Section 4.4. Figure 1 shows the structure for the matrix multiplication data-flow graph, where circles are the tasks and diamonds are the state managers, can be split into the following computational steps:

1. The source (a) of the graph dispatches blocks of matrices A, B, and C to the input nodes.

Listing 1: Simplified MatMul main

The Hedgehog Graph adds these tasks by establishing an edge between them with the method graph.addEdge(SenderNode, ReceiverNode), which will automatically generate FIFO (First In First Out) data queues. Methods from the graph connect tasks to the inputs using graph.input (InputNode) and outputs using graph.output(OutputNode) as shown in Listing 1.

Hedgehog nodes

Hedgehog provides several distinct types of nodes that are used to develop data-flow graphs. The nodes are designed to provide specific functionality in order to achieve a separation of concerns approach; i.e., nodes that specifically target state maintenance, computation, composition, or scalability. Each node type is restricted to consuming multiple input types and producing a single output type, as described in Section 4.3. The node classes form a hierarchy to enable expanding upon the Hedgehog library, as shown in Figure 2 Node is the pure abstract class at the root of the class hierarchy; it is the base class for all Hedgehog node objects. The Graph object represents a computation, it connects a set of consistent nodes and edges to carry out the computation. In Hedgehog, a node must be associated with a parent graph, which binds its constituent nodes to the same device to improve data locality.

The AbstractTask node operates on data; it is abstract as it does not have an implementation for the execute methods of all input types. execute method is called with data that is taken from the task's input queue(s). It executes a computational kernel on the data it receives as shown in Figure 3. This figure shows the execution logic of an AbstractTask, and displays the customization points. The canTerminate function identifies when to stop a task; by default, this is when there are no input nodes connected and nothing in the input queues. This behavior can be modified to break cycles in the graph. The initialize and shutdown functions are called once to carry out pre-or post-computation. Listing 2 for a simplified MatMul task, which implements the execute method to carry out the matrix multiplication operation. The second specialized AbstractTask is the CUDATask. It binds the task's thread to a CUDAenabled GPU, which is done automatically by Hedgehog through a customization of the initialize function. This is used to streamline executing a CUDA kernel, such as by automatically enabling CUDA peer-to-peer access and creating a cudaStream for each instance of the CUDATask.

The last specialized task in Hedgehog is the ExecutionPipeline which is used to simplify multi-GPU programming or other execution contexts. This task creates mirror copies of a graph and distributes data between graph copies according to a developer-specified strategy. Each graph copy is bound to a separate GPU; the graph's data will be local to its tasks and, as such, to the GPU that it binds to. When a StateManager is copied, its State instance is shared across all StateManager instances. This is often used to share data across multiple graphs, and potentially between GPUs. For example, data that is stored within the State can be retrieved and used by another graph from its StateManager. This data may reside in another GPU's memory. It is up to the developer to initiate GPU-to-GPU copying or direct peer access.

Hedgehog is extensible; a specialized task can inherit from AbstractTask to create tasks, for example specialized for AMD GPUs using OpenCL or other libraries.

Nodes input and output data types

Nodes define their input and output types explicitly as template parameters (bound at compiletime) to ensure consistency. The input and output types of two nodes must match in order to connect them by an edge in the Hedgehog graph. Establishing this edge connection will automatically create the queues that are used to store data between the two nodes. A node is restricted to a single output type and can have multiple input types. The input and output types can be anything including types created by a user. For example, if a graph's node requires multiple output types, then these types can be gathered into a single class/struct. The input data types each have their own independent queue. The output of a task can connect to zero or more inputs of other tasks. Data that is sent out of a task is shared with all of its successor tasks. Hedgehog uses smart pointers to avoid the need for a deep copy. Data races that may result from multiple tasks receiving the same pointer are not controlled by the library. It is important for end-users to be aware of their usage patterns on received data.

For example, the accumulate (Acc) task in the matrix multiplication data-flow graph from Figure 1 contains a single input and two output edges. Instantiating this task using Hedgehog is done by specifying P air P C as the input type and Blk C as the output type for the node. An edge is then added from State P C to the Acc, and from Acc to both State C and State P C . The edges will create queues for each of these nodes, if they were not already created, and any data that is produced from the Acc task will be shared with the State C and State P C tasks. In addition, the canTerminate function will need to be customized for the State P C or Acc nodes to break the cycle in the graph as will be discussed in the next section.

Threading model

A graph begins processing data when the graph is executed, which creates a thread for each task in the graph. If a task requests multiple threads (such as MatMul which accepts the number of threads at construction as shown in Listing 2), then custom copies of the task are made using the user-overloaded copy method. The original task and its copies form a task group, where each task in the group is statically bound to a different thread. They share the same input and output edges consisting of queues and synchronization contexts. This forms contention on the input queues, but also offers higher throughput when sufficient data is waiting to be processed.

Hedgehog uses monitor synchronization to operate the edges; a consumer task holds a mutex and will enter the wait state if its input queue is empty. The monitor is implemented with a std::condition variable. It will suspend the execution of the thread if no input data is available. The thread will wake up when notified by another node, this occurs when termination is requested or input data becomes available. When the thread is waiting for data, it will not consume CPU resources. If data is available, then the thread will remove it from its queue. A producer task sends data and signals to its successor consumer task, which transitions the successor consumer task from the wait state to the running state. The running and wait states of a consumer task are managed entirely by the operating system. After a context switch, the signaled consumer task will retrieve the data from its input queue and call its user-defined execute function on the received data. A consumer task can avoid the wait state and potentially the context switch if data were already available in its queue. A consumer task in Hedgehog infinitely fetches data from its input or waits until its canTerminate function returns true; by default this happens when the input connection is terminated and the queue is empty. When termination occurs, the terminated task will signal to all other tasks in the same group to wake-up and terminate properly.

Selecting the degree of parallelism for a task helps improve CPU utilization when executing the graph. For example, in matrix multiplication, the MatMul task can operate with the number of threads equal to the physical core count, whereas the accumulation task can operate with a fraction of the number of physical cores because it is less computationally complex. Specifying a good balance between the number of physical cores and the number of copies of a task can help prevent over-subscription of the CPU for compute intensive tasks, assuming there is enough data waiting to be processed in the task's input queue. In our experience, over-subscription rarely occurs due to the nature of the data-flow approach and how tasks will enter a running state only when data is present in a graph with dozens of independent tasks and edges. This behavior may vary depending on the nature of the algorithm. Furthermore, developers find the idea of tuning the number of threads understandable.

Memory management

Memory management is often needed when GPU computations or limitations due to hardware are involved. Hedgehog implements a memory manager, a tool used by tasks. An instance of a memory manager can be created with a specific data type and then attached to a task. This memory manager instance is thread-safe and limits the amount of data being produced by the task. If the task is copied, such as when building a task group, then all tasks in the group will share the same memory manager.

A memory manager creates a fixed-size pool of user-specified data buffers when the task is initialized. During a task's execute invocation, the task can fetch memory from this pool and will block if the pool is empty. The API provides a mechanism to recycle memory back to the memory manager, which will add memory back into the pool and signal to a task that is waiting. This allows a waiting task to obtain the needed memory and resume execution. The recycle mechanism provides two steps: (1) complete the allocation-deallocation cycle to eliminate memory leaks and (2) update the state of the memory. The state is updated when memory is returned to a memory manager. The state of memory indicates when a memory buffer is ready to be recycled. If the memory is not ready to be recycled, then it will not be added to the memory pool. This functionality is useful to indirectly express memory locality, such as to keep memory resident on GPUs to avoid unnecessary data transfers.

Two types of memory managers exist: static and dynamic. Both managers will fill the pool with objects during task initialization. The dynamic memory manager will call the default constructor for the data objects whereas the static one will call a specialized constructor.

Profiling

Hedgehog provides profiling through a printer class. It is used to represent the current state of the graph. To measure the overhead of profiling, we compare the performance of Hedgehog's first tutorial, the Hadamard product [START_REF]Tutorial 1 -Simple Hadamard Product[END_REF], with and without profiling, and execute it 1000 times on 16k × 16k matrices and 2k × 2k blocks, as shown in Figure 4.

To measure the impact of the profiling, we propose the following statistical analysis.

Let X 1 and X 2 be the execution times with and without profiling, respectively. Given the high number of experiments and the central limit theorem, the experimental averages, X 1 and X 2 , follow a normal distribution. Now define the stochastic variable X as

X 1 -X 2 S 2 1 N + S 2 2 N
. X follows a normal distribution of expectation µ 1µ 2 and variance σ 2 = 1. Let S 2 1 and S 2 2 be estimators of σ 2 1 and σ 2 2 . For the null hypothesis to hold, we propose µ 1 = µ 2 ; there is on average no statistical difference between a computation with and without profiling. X becomes a standard normal distribution.

The size of our sample is N = 1000. For this sample, the values of X 1 and X 2 are respectively x 1 = 1742.28 ms and x 2 = 1743.06 ms and the estimations of the variances σ 2 1 and σ 2 2 are respectively s 2 1 = 14.03 ms and s 2 2 = 12.49 ms. With this sample, the value of X is x = -1.32, whose p-value is 0.4066. These results allow us to accept the null hypothesis.

We can therefore conclude that the average execution runtimes, with or without profiling, do not differ significantly. We believe this is because we gather performance metrics at the node level, which is far less intensive than fine-grained profiling approaches, such as measuring all function invocations. The metrics gathered from the nodes are the execution time, waiting time, and queue usage. Additionally, the graph gathers the creation and total execution times. This information can be presented in various ways depending on options chosen by developers; for example, a task group can be expanded to show the performance of each thread. It is also possible to bind the graph to one or multiple POSIX signals with a graph signal handler to capture them and generate the state of the graph at that instant, such as when a segmentation fault occurs.

A DOT file printer has been developed to generate Graphviz [START_REF] Ellson | Graphviz -open source graph drawing tools[END_REF] DOT files. The developer can use this visualization to understand how tasks and graphs interact with each other and using this representation can immediately recognize the critical path of the computation. Example DOT files can be found at the end of each of the Hedgehog tutorials [START_REF]Tutorial 1 -Simple Hadamard Product[END_REF].

Type checking

Hedgehog uses template metaprogramming to check the data-flow graph consistency at compiletime. It achieves static checking via constructs that use C ++ 17 metafunctions [START_REF] Vandevoorde | C++ Templates: The Complete Guide[END_REF]. This static checking will become simpler and more expressive with C ++ 20 constraints and concepts.

Consistency of smart pointers

Hedgehog manipulates nodes by means of smart pointers based on an RAII (Resource Acquisition Is Initialization) technique, to avoid memory errors in the library. Basically, a smart pointer is an object wrapping a pointer of a given type. In C ++ , smart pointers are modeled by a generic class, shared_ptr<T>, where T is the type of the wrapped pointer.

To accept or reject a node, Hedgehog extracts the node's internal type from its smart pointer and checks this internal type's inheritance.

In Hedgehog, metafunctions have been designed to test the compatibility of two smart pointers.

Node compatibility

Hedgehog checks compatibility rules when two nodes are linked by an edge in the graph:

Output Node 1. If a node is set as the output of a graph, the node and the graph need to share the same output type, as shown in Figure 5a. 2. If a node is set as an input of a graph, the node and the graph need to share at least one input type, as shown in Figure 5b. 3. The output type of the source node must be one of the input types of the target node, as shown in Figure 5c.

I 1 I 1 Output Node X I 1 I 2 (
Listing 3 shows the compatibility test (3), using a static_assert to perform the test and produce a clear error message at compile-time.

static_assert( traits::Contains_v<Output, Inputs>, "The given io cannot be linked together" );

Listing 3: Metafunction invocation testing nodes compatibility

Results

The experiments conducted on Hedgehog help understand the behavior of the library, and the impact that having an explicit representation has on performance.

In our experiments, there are two crucial parameters, the amount of data streaming through the graph (often determined by decomposition parameters such as the block size in a matrix) and the number of threads for each Hedgehog node. The amount of data streaming through the graph affects the degree of parallelism that can be achieved. Because Hedgehog targets coarse-grain parallelism, it is important to determine an optimal decomposition size. For example, if the block size is "too small", then the underlying hardware may have poor utilization because each task does not execute enough instructions compared to the system latency. If it is "too big", then this will reduce the parallelism in Hedgehog because fewer pieces of data feed the tasks. The number of threads for each task will also determine how many elements a task can process in parallel. This can be detrimental if the processor gets oversubscribed. The best methodology for approaching computation in Hedgehog is setting good enough values for a first run. This run is used to generate the DOT file feedback, which identifies bottlenecks, and helps determine better parameters to improve the performance for a specific architecture. These graphs are portable, and only require parameter tuning for new architectures. This methodology justifies experimentation for performance, and was used to help identify the optimal parameters within our results.

Data transfer latency

One of the major costs in our approach is the data transfer latency between tasks. Two parameters affect this latency: (1) the number of threads per task, because all threads in a task group will share the same protected queues, which introduces access contention, and (2) the number of input types, because a task will have one queue per input type.

Figure 6 shows the results of measuring the time to transfer one element between two tasks, averaged over one million data transfers, on a computer with two Intel Xeon Silver 4114 CPU @ 2.20 GHz processors (10 physical cores, 20 logical cores) and 192 GiB of memory. This latency varies between ≈ 1 µs to 10 µs. It grows with the number of threads in the group of receiving tasks and the number of input types. These measurements are low enough that the data transfer costs can be easily hidden by concurrently executing compute kernels. Furthermore, this latency is below that of HTGS for most common cases.

LU decomposition with partial pivoting

One of our experiments was a CPU-based implementation of LU decomposition with partial pivoting [START_REF] Gerson | Study of Exploiting Coarse-Grained Parallelism in Block-Oriented Numerical Linear Algebra Routines[END_REF].

The kernels used in our tasks are mainly composed of calls to the OpenBLAS [START_REF] Xianyi | OpenBLAS: An optimized BLAS library[END_REF] library in single-threaded mode. Our baseline was the LAPACK dgetrf routine, compiled with OpenBLAS in multi-threaded mode, used as a one-off call.

This methodology made it easier to reason about the complex dependencies and to instrument the code in a more sophisticated way that allows the streaming of data in and out of the operation. The stream-based design provides partial fully computed results before the computation completes. This is used to initiate the next step of a computation sooner.

To study end-to-end runtime performance, we ran 10 trials for the Hedgehog and LAPACK (Linear Algebra PACKage) implementations of LU decomposition with partial pivoting on two Intel Xeon E5-2680 v4 CPUs @2.40 GHz (28 physical cores, 56 logical cores) with 512 GiB of DDR4 memory. Figure 7 shows that Hedgehog obtains comparable overall performance to the baseline. In this case, the optimal parameter happens to be a block size of 1024. Additionally, the streaming aspect allows us to start getting partial results over 40 times faster than waiting for the entire matrix to finish computing. This has potential for future research in understanding the intricacies of chaining streaming operations together by composing graphs especially when it can be automated.

Matrix multiplication

Hedgehog is a library that allows computation in a heterogeneous node. We developed a BLASlike GEneral Matrix Multiplication routine (GEMM) and ran it on a node with two Intel Xeon Silver 4216 CPUs @2.1 GHz (16 physical cores, 32 logical cores) with 768 GiB DDR4 memory and four Tesla V100-PCIe with 32 GiB HBM2 GPU, using 128k × 128k single-precision matrices. As shown in Figure 8, the Hedgehog implementation is compared to the NVIDIA counterparts from their GPU-accelerated libraries for basic linear algebra, cuBLASMg and cuBLAS-XT [START_REF]Basic Linear Algebra on NVIDIA GPUs[END_REF]. With an optimal block-size the Hedgehog implementation achieves > 95% of the theoretical peak across 4 GPUs. This final result required multiple iterations that were done using Hedgehog's experimentation for performance. This methodology starts with a simple CPU-only version which is then augmented to use the GPU. The GPU version bundles the data transfers to and from the GPU and the actual computation into its own graph. This was made possible by the composability of the Hedgehog graphs. The GPU graph is further improved by inserting it into an ExecutionPipeline to scale the computation to multiple GPUs. These graphs could be represented as a library, which would simplify the future development of Hedgehog-based applications.

The CPU and GPU versions reused the original representation of state. The primary focus between the versions was optimizing the transfer of data between the CPU and GPU as well as improving the data locality to avoid unnecessary copies. For example, the memory manager was used to stay within the GPU memory limits, and cudaEvents were applied to memory copy events to overlap data motion, computation, and all of the state management that is done within Hedgehog.

These experiments show that the overhead incurred by our approach is not detrimental to performance. It is also amenable to using libraries that target fine-grained parallelism, such as calls to cuBLAS to maximize performance. Lastly, the explicit representation featured in Hedgehog is manageable. The LU decomposition with partial pivoting implementation was completed in the course of one month by a non-domain specialist without any prior knowledge of the Hedgehog library and the intricacies of parallel programming. This experience has continued to indicate to us that having a high-level programming abstraction that can be reasoned about and that maps to execution is extremely valuable.

Conclusion

In this paper, we have presented Hedgehog, a general-purpose library allowing developers to create a parallel computation on a heterogeneous node. It differs from other approaches as it only uses its explicit data-flow representation to implement an algorithm, and relies on data-pipelining to get performance without a global scheduler. This approach operates entirely based on the presence of data and achieves an overhead of ≈ 1 µs to 10 µs when transferring data between tasks in a graph.

Hedgehog supports a separation of concerns approach by providing several distinct components, such as computation, localized state, memory management, and scaling through execution pipelines. The graph is constructed using these components, and its representation is maintained throughout its execution. The developer can use this representation to understand the processing. This helps developers to reason about complex operations at a higher level of abstraction. This was demonstrated through the use of streaming data in and out of an operation to produce partial results as early as possible. This was achieved by explicitly representing the localized states of the computation throughout the algorithm, which aids with understanding when the results are finalized. In addition, Hedgehog can be easily extended as it relies on class inheritance to create new types of tasks, as we have presented for CUDA-based tasks.

Hedgehog operates effectively with heterogeneous computers as well, which is validated by our achieved GPU utilization. This is feasible by making use of Hedgehog's memory management tools, to express memory locality and avoid unnecessary memory copies and keeps GPUs busy.

Future work

Hedgehog has an explicit and static representation of the data-flow graph. The graph structure is known at compile-time. The latest C ++ standard (C ++ 20) brings more tools to perform complex computation at compile-time, such as additions to the constexpr specifier. A compile-time tool will be developed to analyze a Hedgehog graph structure to identify design mistakes such as data-races or cycles that cause deadlock. Another future work will be to extend Hedgehog to operate in cluster environments. One approach is to build one Hedgehog graph per node with direct interaction with external cluster-based communication tools and/or libraries.

Figure 1 :

 1 Figure 1: Matrix multiplication data-flow graph

Figure 2 :

 2 Figure 2: UML simplified class hierarchy of Hedgehog

2

 2 With profiling X 1

Figure 4 :

 4 Figure 4: Hedgehog Tutorial 1 execution time distribution for 1000 experiments on 16k × 16k matrices and 2k × 2k blocks with and without profiling on a Mac Book Pro Mid 2015

Figure 5 :

 5 Figure 5: Compile-time check using metaprogramming techniques

Figure 6 :

 6 Figure 6: Latency analysis for different task's number of threads

Figure 7 :

 7 Figure 7: LU with partial pivoting performance for a matrix of 32768 × 32768 elements

Figure 8 :

 8 Figure 8: Matrix multiplication on a matrix of 64k × 64k elements decomposed in 8k × 8k block elements

  The M atM ul (c) task accepts pairs of A and B blocks, applies matrix multiplication on them, and produces partial results P . 4. The State P C (d) produces a P air P C output when it receives one Blk C and one Blk P that need to be accumulated. 5. The Acc (e) task accumulates the appropriate C blocks with the P blocks. 6. The State C (f) produces an output when Blk C is fully accumulated.

	using MatrixType = float;// Type of Matrix elements
	Order Ord = Order::Column;// Matrix orientation
	// Tasks
	auto matMulTask =
	std::make_shared<MatMul<MatrixType, Ord>>(numberThreadProduct, p);
	auto accTask =
	std::make_shared<AccTask<MatrixType, Ord>>(numberThreadAddition);
	//[...]
	// Build the graph
	matMulGraph.input(stateAB);
	matMulGraph.input(statePC);
	matMulGraph.addEdge(stateAB,matMulTask);
	matMulGraph.addEdge(matMulTask,statePC);
	matMulGraph.addEdge(StatePC,accTask);
	matMulGraph.addEdge(accTask,statePC);
	matMulGraph.addEdge(accTask,stateC);
	matMulGraph.output(stateC);

2. The State AB (b) accepts blocks of A and B, groups them into compatible pairs of blocks, and emits pairs of compatible blocks. 3.

  . The available Hedgehog node types are (1) Graph, (2) AbstractTask, (3) StateManager, (4) CUDATask, and (5) ExecutionPipeline.

	Node
	Receives different types
	of input data, sends one
	type of output data.
	AbstractTask	Graph
	Processes received data.	Organized struc-
	Can be multi-threaded.	ture of nodes.
	StateManager
	Interacts with locked State.
	Restricted to one thread.
	CUDATask
	Specialized Task
	for CUDA GPU,
	bound to a GPU.
	State	ExecutionPipeline
	Shareable thread-safe object.	Duplicates a Graph.
	Useful to express local state.	Customizes routing
		data to the graphs.

  AbstractTask is the StateManager. The StateManager holds a State instance. The State expresses a localized state between two or more nodes. Once an input data is available, the StateManager locks the State it owns and transfers the data to the State. Once the State computation is done, its output data are gathered by the StateManager, the State is unlocked and the gathered data are sent to its connected nodes. The shareable State object contains a mutual exclusion context to ensure it is only accessed by one thread at a time. Both the StateManager and State are fully customizable, but operate with the distinct requirement that StateManager is the only type of task in the graph that should express the state of a computation.

	class MatMul : The first specialized
	public AbstractTask<BlkP, pair<BlkA, BlkB>>{
	private:
	size_t count_ = 0;
	public:
	MatMul(size_t nbThreads, size_t count)
	: AbstractTask<BlkP, pair<BlkA, BlkB>(
	"Product Task", nbThreads), count_(count){}
	void execute(pair<BlkA,BlkB> ptr){
	auto matA = ptr.first;
	auto matB = ptr.second;
	BlkP res {};
	// res <-matA*matB
	cblas_sgemm(matA, matB, res);
	this->addResult(res);
	}
	};
		Listing 2: Simplified version of MatMul task
	start	initialize
	can terminate?	[False]	wait	execute (data)
		notification	[False]
	[True]			unlock
				can	queue
			terminate?
		shutdown		or
				data
			available
					get data
		notify	[True]	[True]	from queue
		all
	terminated		can	[False]	lock
			terminate?	queue
		Conditional			Step
		Customizable conditional	Customizable step
		Figure 3: Task's main runtime steps

Disclaimer

No approval or endorsement of any commercial product by the National Institute of Standards and Technology is intended or implied. Certain commercial software, products, and systems are identified in this report to facilitate better understanding. Such identification does not imply recommendations or endorsement by NIST, nor does it imply that the software and products identified are necessarily the best available for the purpose.


C Farhatt

M Geradint

ON A COMPONENT MODE SYNTHESIS METHOD AND ITS APPLICATION TO INCOMPATIBLE SUBSTRUCTURES

Component Mode Synthesis (CMS) methods are substructuring techniques frequently em ployed in structural dynamics. A given structure is subdivided into components or substructures which are analyzed independently for natural frequencies and for mode shapes. The substructure mode shapes are then assembled to give frequencies and mode shapes of the original structure. In this paper, we construct a substructure interface impedance operator and present a spectral analysis that demonstrates that the method of Craig and Bampton (CB) is the most 'natural' CMS method. Next, we consider the CB method for assembling heterogeneous substructures and recast it into a hybrid variational formulation. We develop fi nite element procedures for 'gluing' non-conforming and incompatible fi nite element substructure models, and discuss their computational aspects. The result is a Hybrid Craig-Bampton (HCB) method that is a fi nite element refi nement of the 'intermediate structure' concept introduced by Hale and Meirovitch, and which can be used as an interface reduction method. It is illustrated with the eigen analysis of heterogeneous and homogeneous fi nite element models of a High Speed Civil Transport (HSCT) aircraft.

t. INTRODUCTION

With the advent of current supercomputers, finite element static analyses of structural problems with 100,000 degrees of freedom (d.o.f.) or more have become almost routine. However, problems with only 50,000 d.o.f. can still challenge the finite element dynamic analyses of structural problems, especially when natural frequencies and mode shapes must be computed. The major difficulty is the expense of computing the eigenpairs of the pencil (K, M), where K and M denote the stiffness and mass matrices of the structure, respectively. Component Mode Synthesis (CMS) is a dynamic substructuring method intro duced by Hurty in 1965 [l], which can be used (a) to reduce the number of d.o.f. in a given fi nite element model in order to reduce the computational cost of a dynamic analysis, or (b) to predict the dynamic behavior of a structural system using results obtained on separate substructures. Since 1965, numerous CMS methods have been proposed. These generally employ coordinate transformations where the gener alized coordinates are defined by a set of sub structure free vibration modes, and may be classified as fixed-interface methods (Hurty (1], Craig and Bampton [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF], Craig and Hale [START_REF] Craig | Block-Krylov component synthesis method for structural mode reduction[END_REF]), free-interface methods (Rubin [4], MacNeal [5], Klosterman [START_REF] Klosterman | On the experimental determination and use of modal representations of dynamic character istics[END_REF]), or hybrid methods (Craig and Chang [START_REF] Craig | On the use of attachment modes in substructure coupling for dynamic analysis[END_REF], Coppolino [8]), depending on whether the modes used are determined with substructure interfaces fully restrained, free of restraint, or partially restrained. Detailed descrip tions of existing CMS methods with their objectives, benefits and major issues can be found, among others, in the review papers by Benfield et al. [START_REF] Benfield | Modal Synthesis Methods. Symposium on Substructure Testing and Synthesis[END_REF], Hintz [START_REF] Hintz | Analytical methods in component modal synthesis[END_REF], and Craig [II], in the textbooks of Craig [START_REF] Craig | Structural Dynamics-An Introduction to Computer Methods[END_REF], Meirovitch [START_REF] Meirovitch | Computational Methods in Structural Dynamics[END_REF], and Cook et al. [START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF], and in the recent monograph by Craig [START_REF] Craig | Component modeling techniques[END_REF]. This paper focuses on the most popular CMS method, namely, the Craig and Rampton (CB) method. Its three main objectives are:

(I) To present a mathematical justification of the CB method which complements its well-known physi cal interpretation and its success in practice. For this purpose, we construct in Section 2 an interface impedance operator whose spectral analysis clearly demonstrates that the set of fixed interface substruc ture free vibration modes augmented with the con straint or junction modes is the 'natural' choice for generalized coordinates. The spectral analysis of the interface operator also leads to the derivation of a new modal participation factor which is economical and simple to implement.

(2) To fully develop a Hybrid version of the Craig and Rampton method (HCB) that is based on a variational formulation (Section 3.1) and which can be used to assemble analysis results obtained on separate incompatible finite element substructure 1 models. In practice, the substructures are often de signed by different teams of engineers, and their respective finite element models often require differ ent mesh resolutions. In such cases, the fi nite element substructure models are typically non-conforming and/or incompatible. The hybrid methodology presented here is a fi nite element refinement of the 'intermediate structure' concept introduced by Hale and Meirovitch [START_REF] Hale | A general substructure synthesis method for the dynamic simulation of com plex structures[END_REF] (see also Meirovitch and Hale [START_REF] Meirovitch | A general dynamic synthesis for structures with discrete substructures[END_REF][START_REF] Meirovitch | On the substructure synthesis method[END_REF]). Its basic ingredients are finite element 'gluing' procedures for assembling heterogeneous substructure models (Sections 3.2 and 3.3). The com putational issues associated with these procedures are fully analyzed (Section 4). We also show that the HCB method inherently reduces the number of inter face coordinates (Section 5), and therefore is a vari ational alternative to the three methods described by Craig and Chang for reducing junction coordi nates [START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF].

(3) To illustrate the efficiency and practicality of the HCB method with the dynamic analysis of realistic three-dimensional structures such as the High Speed Civil Transport (HSCT) aircraft (Section 6).

AN INTERFACE IMPEDANCE OPERATOR AND ITS SPECTRAL DECOMPOSITION

Let { _@(s l}: :::: f• denote a collection of N, substruc tures whose assembly defi ne a global structure .@. Within each substructure, a distinction is made be tween (a) the interface boundary degrees of freedom (d.o.f.) u�l which are common to several substruc tures, and (b) the internal d.o.f. u!'l which exclusively belong to.@<'>. We define the impedance of the global structure as:

Z(w) = K -w1 M (1)
where w is a natural frequency of .@. If every sub structure g&< s l is only subjected to interface boundary forces F�l, the substructure dynamic equations of equilibrium can be written in the frequency domain as: where (2) For a free vibration problem, the internal d.o.f. can be eliminated as:

(3) which leads to the following reduced system:

In eqn (4) above, Z��(w) represents the impedance of the interface boundary of substructure .@ (s) . Since Zfjl(w) is a real symmetric matrix, its inverse Zfjl -'(w) can be evaluated using a spectral decompo sition as:

n = nl"'l cf> (s) cf> (s)T z < s J -1 ( ) = " n n u (1) L.... (s)2 2 n=l Wn -(J) (5)
where cf>�», w�l, and nf'l are respectively the eigenvec tors, eigenvalues, and the size of Zfjl ( w) = Kfjl -w 2 M fjl. Therefore, {cf>�l}: :::: 7f'' and {w�l}:::::f''

are the fixed-interface substructure free-vibration mode shapes and frequencies. The spectral decompo sition of Kfjl-1 is obtained by setting w = 0 in (5) so that one has:

Subtracting ( 6) from [START_REF]A hybrid method of component mode synthesis[END_REF] and using the orthonormal properties of the substructure modes gives after some algebraic manipulations:

and

n = nl'l cf> (s) cf> (s)T z< si -1 ( )-K <si -1 + 2 '°' " " ii (1) -ii (1) L.... (s)2 ( (s)2 2 ) n= 1 (J)n Wn -(J)
Zfjl -1(w) = Kfjl-i + w 2 Kfjl-1 M fjl Kfjl -i Finally, expanding eqn (4) as: and substituting every occurence of Zfjl-1(w) in ( 8) by either one of its two equivalent expressions given in [START_REF] Craig | On the use of attachment modes in substructure coupling for dynamic analysis[END_REF] leads to the following expression of the reduced substructure impedance: (ll) [START_REF] Craig | Structural Dynamics-An Introduction to Computer Methods[END_REF] eqn [START_REF] Hintz | Analytical methods in component modal synthesis[END_REF] above can be re-written in matrix form as:

-

2 bb bi 13 [ M (s) M (s)])[ [<s) J w M�t M!f> .z,�> . ( )
Clearly, cz,� s > is the set of constrained or junction mode, and z�t(w) is the impedance matrix obtained from the well-known Guyan's reduction [START_REF] Guyan | Reduction of stiffness and mass matrices[END_REF]. On the other hand, eqn (ll) shows that z� t (w) is generated by the sequence of rank one matrix products { 4' �lcz,�JT}�::::t> and therefore contains the normal modes.

The above derivation clearly demonstrates that the dynamic behavior of the substructure �< s > can be exactly described with the combination of the junc tion modes and the fixed-interface free-vibration sub structure modes. Therefore, eqns (l)--( 12) constitute a mathematical justification of the Craig and Bampton method.

Remark

If the internal and interface boundary d.o.f. are interpreted respectively as 'slave' and 'master' d.o.f. eqn [START_REF] Benfield | Modal Synthesis Methods. Symposium on Substructure Testing and Synthesis[END_REF] recovers the well-known result that for the lower frequencies (w4Z�t(w) � z�t(w )), the impedance of the reduced substructure is well ap proximated by Guyan's reduction. For the higher frequencies, the fixed-interface normal modes must be represented in the generalized coordinates.

Equation (l l) can be re-arranged as:

n=n(J) z�t(w) = -I ' n= I (14)
which can also be written as:

n =n(-�) z�t(w) = -I ' n=l (15)
Clearly, the contribution of a normal mode to the reduced substructure impedance is dictated by the Several modal participation factors have been suggested in the past (see, for example, Craig and Chang [START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF]), but this topic is continuously revisited (see, for example, Spanos and Tsuha [START_REF] Spanos | Selection of component modes for flexible multibody simulation[END_REF]). The above MPF, 17�> , possesses two attractive features: (a) it explicitly relates the contribution of a given normal mode to the substructure reduced impedance-and therefore to the substructure dynamic behavior, and (b) it is economical and simple to implement as it requires only one matrix-vector product. It may appear at first sight that in order to be able to use the above MPF, 17� s >, one should first compute several substructure modes then retain only those with the highest 17�>. However, if one uses a state-of-the-art eigenvalue solver such as the Lanczos algorithm as implemented by Parlett and Nour-Omid [START_REF] Parlett | The use of a refined error bound when updating eigenvalues of tridiagonals[END_REF], one is able to exploit 17�> during the process of computing the eigenvalues and therefore one would directly compute those eigenvectors with the highest modal participation factors. may be designed using incompatible finite element models. In order to address this class of problems, we develop in this section a hybrid variant of the Craig and Bampton method (HCB) which can be used to 'glue' heterogeneous substructures with non conforming interfaces. The basic idea is to couple the incompatible sub structures with Lagrange multipliers. This has already been discussed by Craig and Chang [START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF][START_REF] Craig | A review of substructure coupling methods for dynamic analysis[END_REF], and by Flashner [START_REF] Flashner | An orthogonal decomposition approach to modal synthesis[END_REF] for the case of compatible models. In [START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF][START_REF] Craig | A review of substructure coupling methods for dynamic analysis[END_REF][START_REF] Flashner | An orthogonal decomposition approach to modal synthesis[END_REF] an algebraic approach appli cable to conforming interfaces is adopted and the Lagrange multipliers are introduced in the formu lation of the problem only to be subsequently re moved with a matrix transformation that reduces the order of the system. In this paper, the Lagrange multipliers are permanent variables of a mixed for mulation. They are introduced in order to enforce an approximate geometric compatibility between the substructures as in the work of Hale and Meirovitch [START_REF] Hale | A general substructure synthesis method for the dynamic simulation of com plex structures[END_REF][START_REF] Meirovitch | A general dynamic synthesis for structures with discrete substructures[END_REF][START_REF] Meirovitch | On the substructure synthesis method[END_REF]. The main contributions of this work are: (a) the development of automatic finite element procedures for 'gluing' non-conforming and/or incompatible interfaces, (b) the discussion of relevant computational issues, and (c) the demon stration of the resulting HCB method with realistic structural examples.

.1. A hybrid formulation

JO( s)

and augmenting it with the inter-substructure conti nuity constraints:

u ( s l =u< qJ s =N5 , q =N s a < s) = u < q ) on u {n < s > ng < ql }. s = 1 . q = I (18)
In eqns ( 17) and ( 18) above, the superscripts (s) and ( q) refer to � < s > and �<q>, respectively, the double dot indicates a second derivative with respect to time, a < s J and £ ( s) are the stress and strain tensors, p < s > is the mass density, u is the displacement field, f is the forcing field which is assumed to act only at the interface boundary d.o.f., and ). ( s , q ) is a Lagrange multiplier function which represents the interface tractions that maintain equilibrium between g < s l and a neighboring g < q> . Since ( 17) is a weak form of the equations of dynamic equilibrium of the substruc tures, we replace the strong inter-substructure continuity equations ( 18) by the following weaker form:

This corresponds to a hybrid formulation where inter-substructure constraints are removed via La grange multipliers (see, for example, Zienkiewicz and Taylor (25]). The specification of ( 17) and ( 18) for static problems corresponds to a saddle-point vari ational principle whose mathematical and compu tational properties are analyzed in Farhat and Roux [START_REF] Farhat | A method of fi nite element tearing and interconnecting and its parallel solution algorithm[END_REF].

Using a standard Galerkin procedure where the displacement and acceleration fields and the La grange multiplier function are approximated by suit able shape functions as: u<s) =Nu<'>; a < s > =Nii('); ;. <s ,q) = Al_ (S) [START_REF] Guyan | Reduction of stiffness and mass matrices[END_REF] eqns ( 17) and ( 19) are transformed into the following algebraic system: [START_REF] Spanos | Selection of component modes for flexible multibody simulation[END_REF] (17) where f < s > is a vector which stores all of the discrete Lagrange multipliers associated with the Galerkin approximation of A. ( s , q>, ns) is the union of all the interfaces between !?}('l and its neighboring substruc tures, and G�l is a vector which stores all of the 'gluing' forces between !?}('l and its neighboring sub structures. Note that if the discretization of the substructure interfaces are conforming, and if purely discrete rather than Galerkin Lagrange multipliers are introduced at every interface boundary d.o.f., the constraint matrix B < sl becomes the signed identity matrix and the second of eqns ( 21) becomes the inter-substructure point-wise continuity equation.

Projecting the substructure equations of dynamic equilibrium ( 21) onto the Ritz basis given by: [START_REF] Parlett | The use of a refined error bound when updating eigenvalues of tridiagonals[END_REF] where p < sl is the vector of normal modes intensities, leads to the following reduced and constrained sys tem: . F (24) and where w is the frequency of the global structure as predicted by the HCB method. The apparent size of the above hybrid system ( 24) is equal to the sum of the number of interface modes, the number of modes retained in every substructure, and the number of Lagrange multipliers that are needed to enforce compatibility between the heterogeneous substruc tures. However, the effective size of ( 24) is much smaller as demonstrated below. Let Z < ;l(w) and Z < ;l' (w) be defined as follows:

0 J [u < •l] n < s ) 2 p � l s =Ns L B < sl u�l = O (
Z < ;6"(w) = Z \;'�(w ) -ro4M b i[n < s > ' -w2 n -1 � r. ( 25
)
Using the above notation, the coupled eigenvalue problem ( 24) can be un-coupled as:

which clearly demonstrates that the effective size of the reduced hybrid system is indeed given by the much smaller size of the operator

[l:�:;f•B < s>z�l'-'(w)B < s> T]-that is, the total number of Lagrange multipliers that are introduced at the substructure interfaces. Further computational aspects related to this operator are discussed in Section 4.

Next, we focus on the approximation of the La grange multiplier function A. < s . q ) via the shape func tions A. First, we consider the case of smooth interfaces where a few Lagrange multipliers are suffi cient to enforce the compatibility between neighbor ing substructures. Then, we address the case of rough interfaces where a larger number of Lagrange multi pliers are needed to glue the substructures. In both cases, we present derivations for problems where the substructures are two-or three-dimensional, but their interfaces are curvilinear.

Low order smooth interfaces

Let e denote the curvilinear abscissae along a smooth interface ns.q) = n<s>nn<q) between substruc tures !?}<s> and !?}<q>, and n��• •' and n �!.., denote the number of degrees of freedom lying on ns , q>nn<s) and ry•q>nn<q>, respectively (Fig. 2). In the case of non conforming interface discretizations, n�! .• , and n � q !.., are not usually equal, and the finite element nodes at one side of the interface do not in general overlap with those at the other side of the interface.

We approximate each interface traction with a directional degree of freedom d(x, y and z displace ments/rotations) with a polynomial of degree p as follows:

k -p A.�• q> = I A.�e k d = 1, . .. , a k -0 ( 27
)
where a denotes the maximum number of degrees of freedom per node. Each of the J x (p + 1) constants

A.� , k = 0, ... ,p, d = 1, ... , Jis an entry ofi .. <s> which appears in [START_REF] Craig | A review of substructure coupling methods for dynamic analysis[END_REF], and can be interpreted as a discrete interface traction between n<s) and n<q>, acting at a location that does not need to be specified by the analyst, and that does not necessarily coincide with a nodal point on ns.q> . where Ia is the J x J identity matrix. As an example, for elements that have two and only two nodes lying on ry•q>nn<s> , and for linear shape functions N., the submatrices L e >��•q>, n = l, 2, k �p, are given by: where e I and e 2 are the Curvilinear abscissae of the first and second interface nodes of element e. The complete element-level constraint submatrices gener ated by an interface element with I nodes lying on ns , q) are given by:

[(e ) B\s , q) ] (e ) B ( s , q) = :

<•>Bj'•q>
B ( s) = L B<s.k) where n<s>nn< k) # <P k B (q) = L B ( q ,k) where n<q>nn< k) # <P (33)
k and can be shown to be sparse matrices.

So far, we have assumed for simplicity that the polynomial approximations of the interface tractions A.�•q> associated with different directional degrees of freedom d have the same degree p. A mathematical discussion on the selection of the polynomial degree p can be found in Dorr [START_REF] Dorr | On the discretization of interdomain coupling in elliptic boundary-value problems[END_REF] who has used a similar 'gluing' procedure for solving the Laplacian problem via domain decomposition. However, a variational interpretation of the constraint equation [START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF] indi cates that p depends on the smoothness of the global mode shapes in the neighborhood of the interface ns , q>. Therefore, p depends on the mode shape to be captured, and, for a given mode shape, on the smoothness of every one of its J field components. In particular, higher order modes may require higher order polynomial approximations of A.�•q>. To this effect, we point out an interesting feature of the procedure outlined above for building the constraint matrices associated with the hybrid formulation of the dynamic equations of equilibrium. Equations [START_REF] Farhat | Using a reduced number of Lagrange multipliers for assembling parallel incomplete fi eld fi nite element approximations[END_REF] show that increasing the degree p of the polynomial approximations of the interface tractions involves only adding a few columns to the existing element level matrices (e)B�•q> and ( e)B�q.s>. Hence, (e)B�•q> and ( e l B�q . s l should be viewed as dynamic matrices which can be efficiently updated, as needed, during the solution of the RHEVP [START_REF] Flashner | An orthogonal decomposition approach to modal synthesis[END_REF]. However, it should be noted that p must always satisfy the condition:

minimum(n��-•l• n�!.. ) ) P < ---�a�--- (34)
which avoids the stiffening of the finest of the two meshes separated by the interface ns.q).

If the desired global mode shapes rapidly vary along the interface ry•q>, a high order polynomial approximation of the interface tractions is necessary because a relatively large number of discrete La grange multipliers are needed to glue the neighboring substructures. However, beyond a certain value of p, the matrices o<s> and B(q l typically become ill-con ditioned, which degrades the performance and accu racy of the numerical solution of the RHEVP [START_REF] Flashner | An orthogonal decomposition approach to modal synthesis[END_REF]. Next, we present an alternative approximation of the interface tractions A.�•q> that is suitable for the case when a rather large number of discrete Lagrange multipliers are needed to interconnect the hetero geneous subsystems.

High order rough interfaces

Here we present a piece-wise polynomial approxi mation of the interface tractions that can accomodate as many discrete Lagrange multipliers as needed to interconnect neighboring substructures, while keep ing the degree of the polynomial approximation unchanged.

First, we estimate a priori the number J x N)., k = 1, ... , N)., d =I, ... , Jwhich form 1<s>. Next, we partition the interface ry.q> into a series of intervals (Fig. 3):

(35) (37)

The first set of equations in (37) imply that A.�q > (ek+.) = A.�t> + i > (ek + 1 ), which shows that A. (s,q> is continuously approximated on ry.q>. The second set of equations in (37) involve the unknown derivatives of the discrete Lagrange multipliers, which have not yet been introduced. We replace these derivatives by their following approximations:

where Aek and A2 ek are defined as:

(38) (39)
Note that (38) requires two additional points, e _, and eN;., Which We Choose as: q> is ap proximated with polynomials, the case where A.� • ql is approximated with piece-wise polynomials requires the analyst to specify the location of the correspond ing physical interface tractions. Therefore from a practical viewpoint, the first approach seems more attractive. However, specifying where an interface traction is needed can be advantageous. For example, if the stress field along I'1 can be qualitatively predicted prior to the analysis, the partition J� • q> can be refined in the areas of oscillation or high concentration, and coarsened otherwise. This would improve the efficiency of the approximation. An adaptive algorithm for locating the intervals J� • q> is given in Farhat and Geradin [START_REF] Farhat | Using a reduced number of Lagrange multipliers for assembling parallel incomplete fi eld fi nite element approximations[END_REF] where the general problem of assembling incomplete finite element ap proximations with Lagrange multipliers is treated.

S= I

Next, we consider the solution of the algebraic system:

(47

)
where uhk+' is the unknown, and ah • is supposed to be known. Such a system arises at every step k of the Subspace Iteration and Lanczos algorithms (Par lett [29]) when applied to the solution of the RHEVP [START_REF] Flashner | An orthogonal decomposition approach to modal synthesis[END_REF]. Expanding the last block of the above equations using the expressions given in (24) leads to:

s =Ns L Bl'> u�Jk+' = 0 (48) S= l which, according to (46), guarantees that -h • + 1TK-h -h •+ I 0 Th ,. h h Vh . . u u > .
ereiore, even t oug A--1s mdefinite, any eigenvalue solver such as the Subspace Iteration or Lanczos algorithms can be still applied to the solution of the RHEVP [START_REF] Flashner | An orthogonal decomposition approach to modal synthesis[END_REF], provided that the starting vectors satisfy the compatibility condition (48).

The solution of the system of linear equations (47)

deserves special attention, as it is the most computa tionally intensive part of any eigenvalue solver. In particular, it would be nice to derive a computational algorithm for solving (47) that requires only substruc ture-by-substructure computations. The hybrid sys tem (47) can be separated into a trivial problem:

(49)

and a more challenging one:

s =N. f L B < s > u�> ' +' = 0. ( 50
)
s= I

At this point, we recall that:

(51) K<i;'f = K�t' if £&<' l is not a floating substructure K<i;' f = K<i:t if £&< s ) is a floating substructure.

The size of the above Lagrangian system is equal to the sum of the total number of Lagrange multi pliers, J x N,, that are introduced to enforce the compatibility between neighboring substructures, and the number of floating substructures, Nr. As noted in Section 3.1, this system is much smaller than the reduced system associated with the interface modes and normal modes intensities. This is why, system (54) above is called here the super-reduced system.

For problems with smooth substructure interfaces, J x N, can be kept small enough so that F, may be explicitly assembled and eqns (56) may be solved with a direct method. For large-scale problems and rough substructure interfaces where N, can be large, the combined direct/iterative domain decomposition method developed by Farhat and Roux [START_REF] Farhat | A method of fi nite element tearing and interconnecting and its parallel solution algorithm[END_REF] is rec ommended for solving the Lagrangian system• (56).

After l k + 1 and cx k + 1 are computed, the remaining components of uh•+i are obtained via simple substi tutions.

S. APPLICATION TO PROBLEMS WITH CONFORMING DISCRETE INTERFACES

l. Modeling and computational savings

The HCB method presented in Section 3 is primar ily motivated by the need for assembling substruc tures with non-conforming discrete interfaces. However, this hybrid method is equally applicable to problems with conforming discrete interfaces and offers attractive computational features over the basic CB method.

Usually, a modal analysis of a proposed structural design is either preceded or followed by a static analysis of that same design. Since mesh generation is by far the most time consuming phase of any finite element simulation (man hours), the same mesh is often used in both analyses. However, experience shows that a finite element mesh tailored for a static analysis is usually much too refined for a modal analysis. In this sense, the weak form of the inter-sub structure continuity conditions ( 19) is desirable even in the case of conforming interfaces. Indeed, the weak form of continuity ( 19) is computationally attractive, as demonstrated below.

For simplicity, consider the case of two substruc tures £&<1) and £&<2l whose finite element discretiza tions perfectly match along a unique interface I'1• The interface equation of the reduced eigenvalue problem resulting from the CB method can be written as: where u b = ul, 1 l = u� l . The interface equations of the reduced hybrid eigenvalue problem resulting from the HCB method can be written as:

Kl, � -Kl,\l Kh1l-1 K�t At each iteration of a given eigenvalue solver, the system of linear equations resulting from (55) can be represented by the operator:

s =N. f A ce -"" K (s) -K (•) K(s)-I K(s) T (57) -f}f~ bb bi ii bi s = l
where the superscript CB stands for the basic CB method, and the system of linear equations resulting from (56) can be represented by the operator:

s =N. f AH CB = L B< s) [K��-K��K!f ) -IK�t]-1BW (58) S= 1
where the superscript HCB stands for the hybrid formulation of the Craig and Hampton method. If discrete Lagrange multipliers are introduced at every 0 0 p < 2)

).

u� l

l p (l) u�> (56) p (2)
).

the inverse of a component of A ce . For large-scale interface systems, a semi-iterative solver such as a Conjugate Gradient (CG) is usually used to solve the corresponding system of equations. Farhat and Roux [START_REF] Farhat | A method of fi nite element tearing and interconnecting and its parallel solution algorithm[END_REF] and Roux [START_REF] Roux | Dual and spectral properties of Schur and saddle point domain decomposition methods[END_REF] have discussed physical and mathematical arguments that establish that CG has a better rate of convergence when operating on a dual operator such as A H ce rather than operating on a primal operator such as A c8 . On the other hand, if polynomial or piece-wise low order polynomial La grange multipliers are introduced on I'1, B<1> and B < 2 > become genuine finite element matrices that have the effect of super-reducing the size of the interface system. For example, if I'1 contains 100 d.o.f. and only 10 Lagrange .multiplier coefficients are needed to glue the two substructures, A ce is 100 x 100, but A" c a is only 10 x 10. The computational advantages of the hybrid formulation are therefore obvious.

Smoothing the mode shapes

Let (v�, <ii �) denote the n th eigenpair of the reduced hybrid eigenvalue problem [START_REF] Flashner | An orthogonal decomposition approach to modal synthesis[END_REF]. If the substructure continuity equations are enforced in a weak form as in [START_REF] Craig | Substructure coupling for dynamic analysis and testing[END_REF], the interface boundary displacement com ponents of the generalized coordinates vector In eqn (59) above, IJI is obtained by assembling the substructure transformation matrices [START_REF] Parlett | The use of a refined error bound when updating eigenvalues of tridiagonals[END_REF]. If the sub structure models have conforming discrete interfaces, it is desirable to compute structural mode shapes that are uniquely defined on every interface nodal point. This can be achieved via the following averaging procedure:

where J!'' is the unsigned counterpart of n< s J_ The matrix product I�:;: f• jj<s> gs)T is a diagonal matrix which computes for each interface degree of freedom the number of subdomains that are connected to it. Therefore, this product should neither be explicitly evaluated, nor inverted.

It should be noted that if the interface boundary components of v�, n = 1, 2, ... , are averaged, the global mode shapes� •• n = 1, 2, ... , may lose their stiffness and mass orthogonality properties. There fore, it is preferable to accept a small discontinuity in the structural mode shapes at the interface bound aries and rather apply the above averaging procedure to the displacement, velocity and acceleration fields when these are evaluated using the computed mode shapes.

NUMERICAL EXAMPLES

Here, we apply the HCB method to the solution of two example problems: (a) an academic eigenvalue problem which illustrates the behavior and basic properties of the proposed hybrid component mode synthesis method, and (b) a modal analysis of a High Speed Civil Transport (HSCT) aircraft which high lights the potential of our method and its important role in re-analysis and re-design cycles.

A two-substructure plate problem is depicted in Fig. 4. Triangular shell elements with 6 d.o.f. per node are used for the spatial discretization. One substructure is clamped at one end. The other one is floating and therefore leads to a singular system of equilibrium. First, we assume that two geometrically non-conforming substructure meshes are given (Fig. 5), and that the first 10 free vibration modes are desired. Three sets of natural frequencies corre sponding to zero (P0), first (P1 ), and second (P 2 )

G-•-•-•- • •ba •----PO +-•-•--P1 a-------P2
order polynomial approximations of the interface tractions l (Lagrange multipliers) are computed. For each case, 10 substructure normal modes are re tained.

In order to check the accuracy of the HCB method, a third 'global' mesh with the same resolution as each of the two substructure models is constructed for the above plate problem, and the first 10 modes are computed via the Subspace Iteration algorithm ap plied to the global problem. The latter solution is here referred to as the 'exact' solution.

All computed results are summarized in Fig. 6. Clearly, if only the first two or three modes are sought after, the P0 and P1 approximations deliver quite accurate results. However, as one would have ex pected, these low order approximations of the inter face tractions do not capture accurately the higher order modes of the structure. On the other hand, the P2 approximation is shown to deliver 10 global modes which are in excellent agreement with the 'exact' solution.

Next, we assume that two geometrically conform ing substructure models are given (Fig. 7). The results reported in Fig. 7 show that the HCB method is equally beneficial to problems with compatible sub structures, since it predicts accurate frequencies at a fraction of the cost of the CB algorithm. Specifically for this analysis, the P2 approximation generates an interface problem of size equal to 18, whereas the CB method would formulate an interface problem of size equal to 66, since there are 11 nodes on the discrete interface between the two substructures. The com puted frequencies are shown to be lower bounds of the actual frequencies as it is proved in [START_REF] Hale | A general substructure synthesis method for the dynamic simulation of com plex structures[END_REF] for this class of problems.

Next, we consider the dynamic analysis of a High Speed Civil Transport (HSCT) vehicle. A structural design and a finite element model for static analysis are initially provided to us. The model includes three (physical) substructures with conforming interfaces and which corresponds to the stiffeners, the skin, and the fuselage of the aircraft. The assembled finite element model is too refined for a modal analysis: it contains 3532 nodes, 8592 triangular shell elements, and 21192 d.o.f. (Fig. 8). Following the arguments presented in Section 5.1, we still propose to use this finite element model for modal analysis but we first partition the global mesh into 4, 8, 16, and 32 submeshes in order to exploit the computational savings offered by the HCB method (Fig. 9). All mesh partitions are obtained with the domain decomposer of Farhat [START_REF] Farhat | A simple and efficient automatic FEM domain decomposer[END_REF].

The results from a static analysis of the above structure under an aerodynamic load have suggested to us that the wing-fuselage junction is poorly con ceived. The lowest natural frequency of the initial structural model was computed and found to be equal to 0.1 Hz, which clearly indicates that the initial design is too flexible. In order to stiffen the vehicle, we have added a carry-through substructure between the wings and the fuselage (Fig. 10). Here, we propose to validate this design modification and illustrate the features of the HCB method with the evaluation of the first six nonzero modes using two different ap proaches. In the first one, the mesh of the carry through substructure is designed to be incompatible and non-conforming at its boundaries with the other parts of the aircraft in order to simulate what really happens when different substructures are designed by different engineers. For this case, the modes of the assembled structure are computed via the HCB method and using a mixture of polynomial and piece-wise cubic spline approximations of the La grange multipliers. Within each substructure, the seven normal modes with the highest participation factor ri�s) are retained. In order to check the accuracy of the 'gluing' procedure, a second mesh that is compatible and conforming at its interfaces with the remainder of the structure is designed for the carry through substructure and the global modes are com puted via the Subspace Iteration algorithm. The obtained natural frequencies are reported in Table 1 for both methods.

Clearly, the 'gluing' procedure is shown to repro duce accurate results. The computational features of the HCB method are highlighted in Table 2 which reports the CPU timings of this dynamic substructur ing method as a percentage of the CPU timings of the Subspace Iteration algorithm. In general as the num ber of substructures increases, the cost for solving the local eigenvalue problem decreases. Unfortunately, this computational gain that is induced by substruc- turing is often offset by the increasing number of interface d.o.f. However, because the HCB method is a genuine interface reduction method, it can be feasible for an increasing number of substructures as reported in Table 2.

Finally, we note that because the finite element mesh corresponding to the original design was given to us, and not constructed by us, we have spent much more effort in designing a geometrically conforming mesh for the carry-through substructure, than we did for constructing the geometrically non-conforming one. In this sense, the HCB method can play an important role in the re-design and re-analysis of structures.

CLOSURE

The Craig and Hampton (CB) method is revisited. First, a mathematical justification of this method is derived via the spectral analysis of an interface impedance operator. It is shown that the set of fixed interface substructure free vibration modes aug mented with the constraint or junction modes is a 'natural' choice for generalized coordinates. The spectral analysis of the interface impedance operator also leads to the derivation of a new modal partici pation factor which is economical and simple to implement. Next, the CB method is recast into a hybrid variational formulation in order to analyze heterogeneous substructures. Finite element pro cedures for 'gluing' non-conforming and incompat ible finite element substructure models are developed, and their computational aspects are discussed. The result is a Hybrid Craig-Hampton (HCB) method that is a finite element refi nement of the 'intermediate structure' concept introduced by Hale and Meirovitch and which can be used as an interface reduction method. The performance of the HCB method is highlighted with the eigen analysis of a High Speed Civil Transport aircraft.

  ) 4' (s) 4' (s)T K (s)T z bb w -L.. w Ww 2 w< s l4 n= I n n If 4' �l is defined as:

  quantity [(K�{-w�>' M �{)4'�>], which physically rep resents the reaction forces induced by the normal mode .z,�> at the interface between �< s > and its neigh boring substructures. Therefore, the above expression (15) leads to a new Modal Participation Factor (MPF): ll(K< '.lw W M < ' >)cz,< s >

3 .

 3 COMPONENT MODE SYNTHESIS VIA LAGRANGE MULTIPLIERS When the substructures { � (s >}�:::: f• are obtained through a partitioning of the fi nite element mesh of the global structure, their finite element models are automatically conforming. However, the beauty of a CMS method resides in the fact that it can be used to analyze a set of substructures that have been designed by different teams of engineers, and to predict the behavior of the corresponding global structure by 'gluing' together the individual substruc ture analyses. In such cases, the submeshes associated with the substructures usually have non-conforming discrete interfaces (Fig. l), mainly because: (a) the corresponding substructures may have different resol ution requirements, (b) the submeshes are often de signed by different analysts, and (c) these submeshes

Fig. I .

 I Fig. I. Substructure models with geometrically non-con forming discrete interfaces.

  Let g < s > , r['>, and n • q> denote, respectively, the volume of the elastodynamic body represented by the substructure �<s>, the portion of its boundary where surface tractions (<s> are prescribed, and the interface boundary between g < s > and a neighboring g < q> . The equations that govern the dynamic un damped response of the global structure can be derived by writing a substructure based virtual work principle:

r

  ( 15 u ( s) p ( s) a ( s) + t5£ ( s) O' ( s» dQ

0 0 B

 0 [START_REF] Craig | A review of substructure coupling methods for dynamic analysis[END_REF] s = Iwhere n < sl ' is a diagonal matrix storing the squares of the substructure frequencies w�l ' , and M < ;l, M < ;? and K �l are the usual quantities obtained with the Craig and Bampton CMS method. Next, considering the free vibration problem and seeking a harmonic sol ution for every field including the Lagrange multi pliers, leads after some algebraic manipulations to the following reduced hybrid eigenvalue problem (Kl,�•) 0 B ( N , ) T 0 0 Q N . , ) 2

Fig. 2 .

 2 Fig. 2. Interface representation.

  and the total contributions of the interface ry.q> to the substructure constraint matrices B<s> and B<q> are obtained via the usual finite element assembly process:B ( s ,q) = L (e ) B ( s , q) e e (I')'• • 'nO«)} B ( q ,s) = L (e )B ( q ,s). ee (I'\'• •) n(l(q)} (32)Finally, if !?}<s> and/or �( qJ have more than one interface with the global structure�. B<s> and B<q> are simply constructed as:

Fig. 3 .

 3 Fig. 3. Interface partitioning.

  and (38) into (37) gives c�k as functions of the unknown discrete Lagrange multi pliers: (41) where </>2k• </>3k• </>4k, '121n '13k• '14k• '2*' '3k• '4k• and V4k are constants that depend only the curvilinear abscissae ekl' ek' ek + l and ek + 2. The B < s > and B < q> constraint operators associated with the piece-wise cubic approximation presented here, are sparse matrices which can be assembled as described in Section 3 . 2. The only difference from the global polynomial case is the building block < •> B� • q> which should be here replaced by its counterpart: Here, the subscript k refers to the interval J� • q> where one edge of the interface element e falls, and � •> f!J� • q> is a K x K submatrix which can be written as: (43) The constants 13l •� i. 13l •> , 13� •� 1, and 13l •� 2 are func tions only of the partitioning of ry . q> and can be computed via the following identification: i A,< s , q> N 'I' -13 < •> 1k-1 M fUk-lAd (rj"•'nrn»ne} + l3l 0> A.�+ 13l:'� 1 A.�+ 1 +13� •� 2A� + 2 • (44) 3.3.1. Remark. Unlike the case where A.� •

4 ,

 4 COMPUTATIONAL ISSUES Now we turn to the numerical solution of the RHEVP (24). Clearly, !(I' is symmetric indefinite and Mil is symmetric semi-definite. Because of the indefin ite nature of !(I', the generalized symmetric eigenvalue problem /(huh= c.0 2 ,.Qhah requires special care. For any given uh' we have: s =NJ s =Ns uh T J(l' ah = L u�) T K bbu�> + L p (s ) T s=l s=l S=Ns x n < sl' p < sl + 2 1.7 L B ( s ) u�). (45) S=I Since K< ;l and Q ( s) 2 are symmetric positive definite, it follows that for any nonzero uh: if L B < '> u�> = 0-+ilhT J(l'ah > O.

  is the reduced stiffness matrix of the substructure � < s > with unconstraint interface boundary d.o.f.-that is, the reduced stiffness matrix of a possibly floating substructure. Consequently, K<i;'� may be singular, and the general solution of the first of equations (50) should be computed as: where K<b' t is a generalized inverse of Kl:�, Rt> stores the basis of the null space of K�� (the reduced rigid body modes), and cx<s)k+' specifies a linear combination of these. The fact that rigid modes do not produce any internal energy can be expressed as: Substituting (52) into (50) for each of the N1 floating substructures leads after some algebraic manipula tions to the following super-reduced Lagrangian sys tem: where s =Ns F, = I s< s )K<i:�• sw S= 1 s=Ns d= I s< s )K<i:U� ) k S= I

Fig. 4 .

 4 Fig. 4. A two-substructure plate problem.

Fig. 5 .

 5 Fig. 5. Discretizations with a non-conforming interface.

0 0interface 0 0

 00 nodal point, B0 l and B (2) become the iden tity matrices, and A ce and A H ce become dual oper ators: basically, every component of A H ce becomes ul, l l 0 b ( J )T p < J)r Kl, 2 � -K�l K ff l-1K�)T 0 b (2)T l ul,

FREOUENCY 0 -Fig. 6 .

 06 Fig.6. HCB results for a non-conforming interface.

0 ��I

 0 v� = ( v � > p <J J • • • v j, N.l p <N, J l]I will be in gen eral discontinuous along I'1• Consequently, the nth mode shape of the global structure computed as: �. = IJlv� where I along the interface boundary r,.

Fig. 7 .

 7 Fig. 7. HCB results for a conforming interface.

FREOIJENCY

Fig. 8 .

 8 Fig. 8. Stiffeners, skin, and fuselage.

Fig. 9 .

 9 Fig. 9. Mesh partitioning (16 substructures).

Fig. 10 .

 10 Fig.10. Carry-through substrucure.

Table I .

 I Eigen analysisRe-designed HSCT structural model HCB Subspace Iteration (16 substructures) Error (%)

	Mode 1	0.572 Hz	0.572Hz	0.0%
	Mode 2	0.690 Hz	0.690 Hz	0.0%
	Mode 3	0.720 Hz	0.720 Hz	0.0%
	Mode 4	1.091 Hz	1.082 Hz	0.8%
	Mode 5	1.450 Hz	1.441 Hz	0.6%
	Mode 6	1.635 Hz	1.612 Hz	1.4%

Table 2 .

 2 Computational costs

	Re-designed HSCT structural model
		CPU (HCB)
	Number of substructures	CPU (Subspace Iteration)
	4	66.6%
	8	35.7%
	16	18.5%
	32	10.8%

Acknowledgements-The first author would like to ac knowledge the partial support of NASA Langley under Grant NAG-1536427, with Dr J. Housner as technical monitor.