Judita Kruse

Clive Thomsen

Rolf Ernst

Thomas Volling

Thomas Spengler

Introducing Flexible Quantity Contracts into Distributed SoC and Embedded System Design Processes *

Increasing design complexity eventually leads to a design process that is distributed over several companies. This is already found in the automotive industry but SoC design appears to move in the same direction. Design processes for complex systems are iterative, but iteration hardly reaches beyond company borders. Iterations require availability of preliminary design data and estimations, but due to cost and liability issues suppliers often hesitate to provide such preliminary data. Moreover, companies are rarely able to judge the accuracy and precision of externally estimated data. So, the systems integrator experiences increased design risk. Particular mechanisms are needed to ensure, that the integrated system will meet the overall requirements even if part of the early estimations are wrong or imprecise. Based on work in supply chain management, we propose an inter-company design process that is based on formal techniques from real-time systems engineering and so called flexible quantity contracts. In this process, formal techniques control design risk and flexible contracts regulate cooperation and cost distribution. The process effectively delays the design freeze point beyond the contract conclusion to enable design iterations. We explain the process and give an example.

Introduction

Efficient design processes and new strategies are needed to meet the challenges of the increasing complexity of current and future embedded system and SoC designs. Industries like automotive engineering and space applications have established sophisticated supplier-integrator chains. We expect that the semiconductor SoC design will follow the same path in SoC development. The success of IP providers such as ARM and their frameworks for intercompany exchange and interoperability of documents and tools [START_REF]SPIRIT Consortium[END_REF] are clear indicators Such supply chains turn out to be the first step towards a distributed design process. Today we already find distributed design processes in automotive engineering and space applications, where software plays an important role. Distributed design processes are spread over several companies and cover the whole design flow from subcomponents to the integration of the overall system.

Increasing SoC complexity, multiple layers of hardware dependent software, application programmer interfaces, and third party software libraries give raise to further design process distribution over a heterogeneous set of players that follow divergent business goals. They must be coordinated via engineering and business mechanisms. It is about time for systematic approaches to that issue that put these mechanisms into context.

Especially in complex product development implementation and integration are done in the late phases of classical design flows (e.g. specified for ESA projects by ECSS standard 'System Engineering' [START_REF]ECSS-E-10A Standard, System engineering[END_REF]). Hence accurate assertions about design data and objective requirements are late as well. But contracting based on design specifications containing requirements is done at the beginning of a new design. Therefore integrator and suppliers estimate their needs conservatively leading to increased cost and suboptimal products.

One key problem in distributed design processes is to reach true collaborative inter-company developments, while keeping corporate know-how proprietary at the same time.

The approach in this paper was motivated by the SpeAC project in which a large space systems company, EADS Astrium, develops an inter-company design flow with two suppliers and applies it to a concrete satellite hardware/software subsystem development as a demonstrator. As a prerequisite, design data are shared between the companies using restricted access mechanisms to protect IP and management internals.

One of the side effects of this inter-company design flow is the availability of formalized design data using standards such as SystemC [START_REF]SystemC[END_REF] for hardware description or SysML [START_REF]Systems Modeling Language[END_REF] for overall system modeling including thermal, power, and mechanical design aspects. It has been shown in the project that the necessary data for the process proposed in the following paper can be formulated and communicated in SysML.

After the discussion of related work in chapter 2, we introduce flexible quantity contracts and their application to distributed SoC and embedded system design processes in chapter 3. We will illustrate our approach with an example given in chapter 4 and conclude in 5.

Related work

The realization of a continuous SoC and embedded system design flow is of great concern and is accelerated by the corporate as well as by the science community. Consortia like the VSI Allicane [START_REF]Virtual Socket Interface Alliance[END_REF] and SPIRIT [START_REF]SPIRIT Consortium[END_REF] are founded to enable the development of SoCs with a special focus on configurable predesigned IP-blocks. They specify a catalog of standards, listing essential design data for different groups of components.

Engineering science discusses different approaches to enable an efficient SoC design process like platform or component based strategies [START_REF] Chang | Surviving the SOC Revolution[END_REF][START_REF] Keutzer | System-level design: Orthogonalization of concerns and platform-based design[END_REF][START_REF] Kopetz | Component-based design of large distributed real-time systems[END_REF][START_REF] Zhang | Component selection and matching for IP-based design[END_REF]. Formal realtime analysis techniques are closely connected to these approaches and address heterogeneous systems especially [START_REF] Chakraborty | A general framework for analysing system properties in platform-based embedded system designs[END_REF][START_REF] Pop | Bus access optimization for distributed embedded systems based on schedulability analysis[END_REF]. To enable advanced real-time analysis techniques particular system models are introduced in [START_REF] Richter | A formal approach to MpSoC performance verification[END_REF][START_REF] Palencia | Best-case analysis for improving the worst-case schedulability test for distributed hard real-time systems[END_REF]. In our paper real-time requirements are used as an example of nonfunctional system properties. This class of systems is well suited for illustrating problems emerging from the interference of different components.

Requirement management and tracing is an important issue in design processes of complex products. The commercial tool DOORS [2] is widely used in the system industry to achieve this, while the eurostep AP233 demonstrator [1] deals with the task of mapping requirements to system components. Booth are used in the mentioned development project. The interaction of competing requirements is not addressed in this paper. For an overview of current approaches refer to [START_REF] Dardenne | Goaldirected requirements acquisition[END_REF].

Contractual agreements have been studied in a variety of academic disciplines, notably economics, engineering, law, and operations management. In contrast to various qualitative contributions (e.g. [START_REF] Cooter | Law and economics[END_REF][START_REF] Kawasaki | The design of contracts: evidence from Japanese subcontracting[END_REF][START_REF] Telser | A theory of self-enforcing agreements[END_REF]), quantitative treatment can be found in operations management, particularly in work analyzing supply-chains. Here, settings of at least two independent actors, interlinked by flows of information, physical goods, and financial funds are analyzed [START_REF] Christopher | Logistics and supply chain management[END_REF]. Common premises include asymmetric access to information as well as significant difficulties to monitor other actors' performance (hidden action) [START_REF] Tsay | Modeling Supply Chain Contracts -A Review[END_REF][START_REF] Wang | Contracting Structures for Custom Software Development: The Impacts of Informational Rents and Uncertainty[END_REF].

In this context a central phenomenon causing inefficiencies is titled double marginalization and arises in supply-chains where decisions of individuals are driven by their subjective perception of the supply-chain incentive structure. This does not necessarily correspond to the optimal decision basis as compared to a central decision maker [START_REF] Corbett | Designing supply contracts[END_REF]. A favorable contractual agreement in this context is one, that at the same time promotes efficiency, Pareto optimality, and incentive compatibility, hereby incorporating individual rationality [START_REF] Whang | Contracting for Software Development[END_REF]. A comprehensive review of contract analysis in the field of supply-chain management is provided by [START_REF] Tsay | Modeling Supply Chain Contracts -A Review[END_REF]. First analysis of contractual issues in software business has been proposed by [START_REF] Richmond | Software development outsourcing contract[END_REF][START_REF] Whang | Contracting for Software Development[END_REF].

For situations of uncertain demand the managerial flexibility to change a chosen course of action is usually associated with better performance [18]. The objective of flexible contracting in this sense is to increase the flexibility of certain entities. Accordingly uncertainty and thus risk is being shared with respect to the contractual agreement in place. Of particular interest for the research presented are contracts with inbuilt flexibility in terms of quantity, so called flexible quantity contracts [START_REF] Tsay | The Quantity Flexibility Contract and Supplier-Customer Incentives[END_REF]. Hereby contractual clauses exist, that define conditions under which the exact amount to be purchased may diverge from an initial estimate. Instances of flexible quantity contracts include total minimum and buy-back contracts. However, in order to adopt flexible contracting schemes to the development of microelectronic devices, specific aspects of interdependent design parameters need to be incorporated. To our knowledge yet no contribution exists in this field of analysis.

Introducing flexible quantity contracts

Fixed price versus flexible quantity contracts

Figure 1a outlines a classical distributed design process based on fixed contracts. In case of SoC or embedded system design the system S is concretized by a hardware/software system under development.

Initially the integrator performs system partitioning and optimization. Subsequently she formulates her needs regarding the components s i as requirements r i and releases final component specifications at a certain point in time t * called design freeze. Requirements concern functional and non-functional design data and can originate from the systems environment, from technical needs, or from customer requests. Requirements can be classified according to critical and non-critical requirements, depending on whether the system must or should achieve the associated demands. In fixed contracting scenarios integrators will pass on critical system requirements to components as critical requirements again, to avoid the risk that a critical overall system requirement is missed (see figure 1b).

In a next step the suppliers prepare offers that correspond to the requirements predicting their design data and cost with special emphasis on critical requirements that must be estimated very conservatively. The offer includes these predictions that are passed to the integrator as higher level component assertions a i . With contracting at t 0 the suppliers ultimately give fixed assertions that their implementations will meet the requirements (figure 1c). By fixing the price p i of a component s i , the specific price associated with its assertion a i is determined by

p i = p(s i ) = p(a i )
In case the integrator refines her requirements during contract negotiations, design freeze and contracting coincide at t * = t 0 . Regarding a simple single requirement system consisting of n components and assuming that each component s i is described exactly by one assertion a i , the price P is given at t 0 by

P = n i=1 p(s i ) = n i=1 p(a i )
In a subsequent implementation phase, the suppliers take the risk of not complying with assertions, while the integrator's risk is that of changing requirements due to unexpected contingencies, or even worse, of a faulty integration. In case of the necessity to change requirements, usually formal change requests have to be filed.

With growing complexity a series of drawbacks of fixed price contracts can be identified. As subsequent changes of requirements are not explicitly part of the contract (as opposed to the procedure how changes are dealt with, i.e. change management) they consequently are neither foreseeable in terms of cost nor of performance and development time. In this situation rational behavior dictates for the integrator to build in more safety in the system architecture and thus in request for proposals than initial analysis suggests. In addition suppliers benefit from tight requirements, as these represent extra business from their point of view, but at the same time due to an increased exposure to development risk raise their safety margin. The inevitable consequence is the pursuit of over-dimensioned designs.

In fixed contracts potential design trade-offs at system level are limited to the pre-contracting phase (t < t 0 ), resulting in an early design freeze and consequently hampering reactions promoted by better information available. Also trade-off cost information will be private to the individual subjects and thus only, if at all, available to the other actors after renegotiation. Accordingly when conducting design trade-offs, decisions are restricted to local information. Results therefore solely reflect a fraction of the design project's decision basis. This situation corresponds to the mentioned double marginalization phenomenon in supplychain management. Similar contracting schemes should therefore be applicable.

In the problem setting described, flexible agreements offer an opportunity to both accelerate the decision process and increase its efficiency. The basic idea of flexible contracting in distributed decision making is to on the one hand defer design freeze by extending the period for possible design refinements and on the other to increase the efficiency of decentralized decision making. A basic approach of setting up efficient decentralized process control, as proposed in economics, is to impose the cost structure of the total value-chain (i.e. the design-consortium) on the decisionmaking entity (i.e. each individual actor), thereby internalizing external effects [START_REF] Mathewson | An Economic Theory of Vertical Restraints[END_REF]. One commonly accepted approach in supply-chain management lies in designing more complex contracts by means of contractual baseline and incentive schemes applied [START_REF] Tsay | Modeling Supply Chain Contracts -A Review[END_REF].

As a prerequisite of adopting a contracting scheme of such kind to design processes, initially suitable measures substituting quantities of physical goods in supply-chain applications need to be identified. In the following we will introduce structured estimation data as a capable measure for flexible contracting.

Flexible contracts in design processes

To enable flexible contracts in SoC and embedded system design, we propose a structured estimation format for assertions and introduce set-critical requirements. Assertions in the structured estimation format consist of three values: a conservative estimation a gua , a target a tar , which denotes the baseline, and a best case a top , replacing physical quantities.

The term set-criticality will be used for requirements which are critical in conjunction with other requirements, but are not critical by themselfs. This statement can be held even for requirements derived from critical requirements, provided that several combinations of interdependent components with different implementations exist. Obviously any kind of optimization benefits from this fact. Doing so a critical system requirement can be attenuated to set-critical component requirements.

A distributed design process based on flexible contracts is illustrated by figure 2a. The integrator starts with system partitioning and component specification. But in contrast to fixed price contracts she will pass a critical system requirement as set-critical requirements to derived components (figure 2b). The suppliers submit offers with assertions in the structured estimation format a i = (a gua i , a tar i , a top i ) instead of delivering only a fixed assertion (figure 2c), along with a corresponding set of prices p i = (p gua i , p tar i , p top i ). Using the flexible contracting scheme comprehensively illustrated in figure 3 the integrator will determine the slope corresponding to linear interpolations. In order to compensate the supplier for the increased exposure to development risk, an additional premium b i will be paid for deviations to the initial target. From the integrator's point of view, b i is used as an incentive for the suppliers to direct the development of the system in a suitable way. Applying sensitivity analysis she can identify which components have the largest influence on the critical property and can assign them a higher b i .

Regarding a single requirement and assuming that each component s i is described exactly by one assertion a i , the specific price p i associated with a i is determined by

p i (a i ) = p tar i + (v gua i -b i ) • (a i -a tar i ) for a i < a tar i p tar i + (v top i + b i ) • (a i -a tar i ) for a i ≥ a tar i
Based on the sets of assertions given by the suppliers, the integrator executes a sanity check, to assure that the critical system requirement will not be violated, according to integrator defined rules. Examples for sanity check metrics are:

The critical system requirement r with a safety margin of m% will be achieved (with x% of a top i -a gua i ) if. . . one-sided: . . . each final component implementation will change the target assertion a tar i at most by x%. conservative: . . . one-sided, but at most y final component implementations will not improve the guaranteed assertion a gua i . optimistic: . . . one-sided, but at least y final component implementations will reach the top assertion a top i . balanced: . . . y of n final component implementations will improve the target assertion a tar i at least by x%, while (n -y) final component implementations will fail the target assertion a tar i at most by x% . weighted: . . . the final component implementations will change the target assertion a tar i by an individual coefficient of x%, e .g. depending on assessments based on former experiences with the particular supplier.

If the sanity check was passed successfully, contracting is done at t = t 0 with the offered assertions and their corresponding flexible pricing schemes.

During the subsequent prototyping phase, the suppliers deliver n iteratively refined estimations at t = t i for the final implementation, but subjected to the boundaries given by a gua and a top , and before a certain point in time t * . At design freeze the final specification will be determined, and the so far set-critical requirements become critical ones with

r i = a * i for a price of P * i = n i=1 p(a * i ).
With flexible contracts an efficient strategy is proposed to establish system wide risk management and optimization within the refinement loop of distributed SoC and embedded system design processes. Change requests can be omitted for the most part, because possible variations of implementations have been anticipated by deferring design freeze beyond t 0 , the point contracting is done.

Example

To test the applicability of our approach, we calculated different scenarios for the example system shown in figure 4. An integrator assigned to design the system partitioned it choosing the following architecture: three functions F 1 , F 2 , F 3 should be implemented in software, running on a customized CPU. Due to hard deadlines with periodic activations given for the three functions a realtime scheduler is needed. The critical system requirement is r critical ='CPU-utilization ≤ 100% and can be verified performing schedulability analysis provided, e.g. by SymTA/S [START_REF] Symta | [END_REF].

The software shall be subcontracted to a software supplier, and the implementation of the CPU to a hardware supplier respectively. The suppliers include the estimations of the worst case execution times for the software and the clock rate for the CPU along with prices (in mu, i.e. 'money units') as assertions in their offers. Example data are shown in table 1, with the calculated slopes v gua i and v tar i . The integrator chooses the scheduling policy a sched ='simple static priority preemptive' [START_REF] Liu | Scheduling algorithm for multiprogramming in a hard-real-time environment[END_REF]. Depending on the risk she is willing to take, the integrator will pick an appropriate sanity check at t = t 0 . As the results in table 2 show, the selection is crucial. With an assumed safety margin of 5% the critical requirement 'CPU-utilization' must not exceed 95%. In our example the integrator chooses the balanced check, which passed successful. To find reasonable premiums b i , a sensitivity analysis tool, as offered by SymTA/S is valuably. Results a sens i of the sensitivity analysis and the chosen b i are shown in table 3.

After contracting but before design freeze the suppliers start prototyping and update their estimations frequently for a predefined number of n iterations. As design progresses, the integrator will monitor expected system performance, based on estimations with constantly increasing accuracy. For illustration table 4 shows data for three hypothetical iterations at t 0 ≤ t 1 , t 2 ,t 3 ≤ t * .

When validation stops, the integrator freezes the design. In the example the critical system requirement was met and a fairly good CPU-utilization was reached. Therefor the integrator will tighten the component requirements r i from set-critical to critical ones according to the assertions at t 3 . The associated prices can be found in table 5 with calculated prices for a fixed price scenario. The fixed price values are based on the assumption, that the integrator would have called for bids on critical requirements equal to the target implementation, which leaves a significant safety margin regarding CPU-utilization. The comparison shows that a fair value-chain was obtained with a win-win situation for all actors. Supplier of F 1 receives slightly less, but avoids change requests due to a missed critical requirement. Supplier of F 2 is rewarded for her effort, while supplier of F 3 earns significantly less, though is prevented from potential contractual penalty. From the integrator's perspective the development of the system is less expensive without an oversized dimensioning.

Table 2. Sanity checks

Employing flexible contracts with estimations the integrator gains an extended period for design trade-offs. Consequently cost reductions can be achieved from the integrators perspective. Considering the suppliers, an overall loss in business can be identified in some cases. However incentive compatibility is created, because an extra premium is paid, due to the compensation scheme. Moreover broader requirements are likely to reduce the suppliers internal costs. Consequently applying the flexible contracting scheme an enhanced design process is facilitated, complying with Whang's criteria for optimal contracting [START_REF] Whang | Contracting for Software Development[END_REF].

Conclusion

Inter-company design processes gain importance. Contractual constraints and IP protection issues entail an overly conservative design style. Motivated by a satellite component design project, a design process based on flexible quantity contracting and formal analysis was proposed that reduces overdesign but allows to control design risk. Redesign cycles can be reduced, because critical component requirements can be traded off against each other. A simple example shows the incentives for the different players in the design process and demonstrates fairness.

Figure 1 .

 1 Figure 1. Flow with fixed price contracts

Figure 3 .

 3 Figure 3. Flexible pricing scheme

Figure 4 .

 4 Figure 4. Model of the example system

  check x[%] y pass info one-sided 0y 75.24% CPU-util. one-sided -5y 85.4% CPU-util. one-sided -10n 105.77% CPU-util. conservative 0 1 n failed for CPU optimistic -10 1 y balanced 10 2 y n! y!•(n-y)! analysis steps weighted 1)y 85.4% CPU-util. 1) x(F 1, F 2, F 3, CP U)[%] = (-20, +10, -50, 0)

F1

  aF 1,F 2,F 3[cycles], aCP U [MHz], pi[mu]

Table 1 . Example data

 1

		a gua i	p gua i	a tar i	p tar i	a top i	p top i	v gua i	v top i
	F1	30	10 20 50 10 150	-4 -10
	F2	60	10 40 60 20 180 -2.5	-6
	F3 150	10 100 60 80 120	-1	-3
	CP U	50 100 100 200 200 600	2	4
	aF 1,F 2,F 3[cycles], aCP U [MHz], pi[mu], vi[ mu [a i ] ]

Table 3 . Sensitivity analysis based on a tar

 3

	i

Table 4 . Risk management

 4

	t aF 1 aF 2 aF 3 aCP U CPU-util.[%]
	t1 21 55 120 100	91.95
	t2 28 42 135 100	94.57
	t3 21 37 150 100	88.52
	aF 1,F 2,F 3[cycles], aCP U [MHz]

Table 5 .

 5 Pricing

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

Acknowledgement

We would like to thank Roland Müller and others from Astrium, Hubert Stich from EuroTelematik, and Peter Ganal from Tecnotron for their valuable input and feedback from design practice.

* This work is supported by a grant from EADS and from the EU in the SpeAC project (MEDEA+ A508).


Mario R Casu
email: mario.casu@polito.it

Luca Macchiarulo
email: lucam@hawaii.edu

A New System Design Methodology for Wire Pipelined SoC

Wire Pipelining (WP) has been proposed in order to limit the impact of increasing wire delays. In general, the added pipeline elements alters the system such that architectural changes are needed to preserve functionality. We illustrate a proposal that, while allowing the use of IP blocks without modification, takes advantage of a minimal knowledge of the IP's communication profile to dramatically increase the performances. We showed the formal equivalence between WP and original system and proved the higher performance achievable through a relevant case study.

Theory and Implementation

A signal is a set of events e, i.e. couples e = (v, t) of values v and tags t. In our sequential systems, tags can be thought as clock ticks. Processes exchange signals by means of channels. If we add wire pipeline elements to channels, these initially hold void data that we denote using τ symbol. The behavior is modified such that sequences of valid events, are interrupted by void values. A generic realization of a channel in time interval [t 1 , t N ] is (v 1 , t 1 ), τ, τ, (v 2 , t 2 ), . . . , τ . . . Let us filter out void symbols τ and find the maximum tag N such that every signal has a sequence of at least N values. If all sequences are identical to the original system from 1 to N , the two systems, with and without WP, are N-equivalent. If they are equivalent ∀N , they are said equivalent. A wrapper must enclose the processes in order to guarantee equivalence and perform the following tasks: 1) τ -filtered inputs are buffered with semi-infinite fifos and stored in positions progressively numbered with the tags. 2) a synchronizer keeps trace of the current tag and, a) when all identically tagged inputs are available (synchronous signals) dispatches them to the internal process, removing them from the fifos; b) if at least one input does not have the current tag, the process is stalled.

3) In correspondence with the stall, τ is sent to all outputs. This wrapper could be simplified so that all valid signals ( = τ ) already processed (and so eliminated from the fifo) are counted. The synchronizer's task is to trigger the process if all inputs with tag equal to the counter are present.

If we suppose that some processes are such that, for some of their internal states, not all inputs are simultaneously read, then it is possible to advance the computation even before all the corresponding inputs have reached it (relaxation of synchronicity). This is the key for an overall increase of performance. An analogous observation is in [START_REF] Singh | Generalized Latency Insensitive Systems for Single-Clock and Multi-Clock Architectures[END_REF] without proof of equivalence. We introduce a new element in our wrapper: an oracle that decides which inputs are needed for the next computation. The new wrapper's features are: 1) if all inputs required by the oracle are present, the computation is triggered and the fifo updated.

2) The synchronizer discards all inputs whose tag is smaller than the counter value ("old" tags).

Synchronicity is guaranteed by the fact that the fifo actually discards valid non-necessary inputs which became useless due to the process' blindness to them.

The step to make the entire consideration practical is making the fifo length finite. Thus, back-pressure has to be implemented with a stop signal propagated back to the wire pipeline elements that, in order not to lose data, should contain, together with the pipeline register, at least one auxiliary register for saving a valid incoming datum when a stop is received. A simple FSM already described in [START_REF] Carloni | A Methodology for Correct by Construction Latency Insensitive Design[END_REF] and called relay station (RS) suffices to guarantee correct functionality. When also the auxiliary register is full, the stop is propagated to the previous RS up to the source process.

Instead of counting the absolute tag, it is sufficient to use an initialized counter that records the lag between the most recent and old data received by each channel. From the practical standpoint, it is also not necessary to send the tag together with the signal, but only a bit indicating its validity [START_REF] Carloni | A Methodology for Correct by Construction Latency Insensitive Design[END_REF], as the distributed counters will take care of the tag values, thanks to the ordering properties of the signals.

Wrappers with and without the additional oracle that uses the processing information were described in VHDL and simulated. We evaluated the wrappers' area with several synthesis experiments on a 130 nm technology. The overhead was always less than 1% with respect to an IP of 100 kgates. As for the timing of the logic used within the wrappers, it was never critical in all our experiments.

Case Study

We considered a processor made out of five components to be enclosed within our wrappers and with pipelined connections as shown in fig. 1 pitfalls are the netlist loops. Let us define the throughput Th as the average number of valid data produced in every clock cycle. A loop containing m process and n pipeline delays has Th= m m+n in case wrappers do not implement oracles. The worst loop dominates the system Th. The netlist loops are shown in fig. 1. We built the system with a minimal instruction set and used two basic programs to cover the spectrum of applications: a strictly data dependent problem, extraction sort, and a matrix multiplication. The processor was modelled with a RTL VHDL description and implemented in two different fashions, multicycle and pipelined.

The migration to a standard WP system, WP1 in our experiments, is straightforward: blocks are encapsulated by wrappers and wires segmented by relay stations. In the case of our new technique, WP2, the input signals are complemented, when possible, by a processing signal derived from the process operation. However, in our case (which we think representative of practical cases), the effort was minimal.

Experimental Results

We ran the following experiments for both multiprocessor and pipelined case but for space reasons only the latter case results are reported in table 1: 1. The golden system is simulated. The number of clock periods is used for the computation of the WP's performance. The throughput without WP is of course 1.0. 2. We added single relay stations (RS) for each and every connection in the system's topology (rows 2-11).

3. We built a system with 1 RS in each connection (row 12). 4. For the Matrix Multiply only, we repeated the simulations for at least 1 RS everywhere and 2 somewhere (rows 13-22) and 2 RS everywhere but in some links (rows 23-25).

We would like to draw some overall conclusions: 1. All results are in favor of the proposed WP2 system. 2. The advantage depends on the features of the communication channel at stake: In the Multicycle case, not reported in table, the CU-IC loop is excited only every 5 cycles due to the non pipelined sequence of 5 phases Instruction Fetch, Decode and contextual Operand Fetch, Execution, Memory access and Write-back. That's the reason of the best improvement of WP2 in this loop with respect to WP1 (60%) where such loop is statically present at all time. Other channels accessed more frequently give less advantage. 1. Even if the computations are tighter in the pipelined case (in the same clock cycle different loops are exercised at the same time) we still observe the relevant advantages of WP2.

Figure 1 .

 1 Figure 1. Case study.

Table 1 . Sort: Pipelined Case.

 1

				Th	Th	WP2 vs.
		RS Configuration	Cycles	WP1	WP2	WP1 (%)
			Extraction Sort
	1	All 0 (ideal)	1559	1.0	1.0	0%
	2	Only CU-RF	2078	0.75	0.75	+0%
	3	Only CU-AL	2090	0.667	0.75	+13%
	4	Only CU-DC	2078	0.75	0.75	+0%
	5	Only CU-IC	3118	0.5	0.5	0%
	6	Only RF-ALU	1886	0.667	0.83	+25%
	7	Only RF-DC	1569	0.667	0.99	+49%
	8	Only ALU-CU	1681	0.667	0.93	+40%
	9	Only ALU-RF	1701	0.667	0.92	+38%
	10	Only ALU-DC	1634	0.667	0.96	+44%
	11	Only DC-RF	1624	0.667	0.96	+44%
	12	All 1 (no CU-IC)	2325	0.5	0.67	+34%
	13	Optimal 1 (no CU-IC)	1952	0.667	0.80	+20%
			Matrix Multiply
	1	All 0 (ideal)	2778	1.0	1.0	0%
	2	Only CU-RF	3704	0.75	0.75	+0%
	3	Only CU-AL	3724	0.667	0.75	+13%
	4	Only CU-DC	3704	0.75	0.75	+0%
	5	Only CU-IC	5557	0.5	0.5	0%
	6	Only RF-ALU	3596	0.667	0.77	+16%
	7	Only RF-DC	2824	0.667	0.98	+47%
	8	Only ALU-CU	2851	0.667	0.97	+46%
	9	Only ALU-RF	3436	0.667	0.81	+22%
	10	Only ALU-DC	3058	0.667	0.91	+37%
	11	Only DC-RF	2994	0.667	0.93	+40%
	12	All 1 (no CU-IC)	4703	0.5	0.59	+18%
	13	All 1 and 2 CU-RF	4775	0.5	0.58	+16%
	14	All 1 and 2 CU-AL	4703	0.4	0.59	+48%
	15	All 1 and 2 CU-DC	4703	0.5	0.59	+18%
	16	All 1 and 2 CU-IC	8335	0.33	0.33	0%
	17	All 1 and 2 RF-ALU	5521	0.4	0.50	+25%
	18	All 1 and 2 RF-DC	4703	0.4	0.59	+48%
	19	All 1 and 2 ALU-CU	4776	0.4	0.58	+45%
	20	All 1 and 2 ALU-RF	5235	0.4	0.53	+33%
	21	All 1 and 2 ALU-DC	4919	0.4	0.56	+40%
	22	All 1 and 2 DC-RF	4919	0.4	0.56	+40%
	23	Optimal 2 (no CU-IC)	4919	0.4	0.56	+40%
	24	All 2 (no CU-IC)	6555	0.33	0.42	+26%
	25	All 2 and 1 CU-RF	6555	0.33	0.42	+26%

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE


Minas Dasygenis

Erik Brockmeyer

Bart Durinck
email: bart.durinck@imec.be

Francky Catthoor
email: francky.catthoor@imec.bedimitriossoudrisdsoudris@ee.duth.gr

A Memory Hierarchical Layer Assigning and Prefetching Technique to Overcome the Memory Performance/Energy Bottleneck *

The memory subsystem has always been a bottleneck in performance as well as significant power contributor in memory intensive applications. Many researchers have presented multi-layered memory hierarchies as a means to design energy and performance efficient systems. However, most of the previous work do not explore trade-offs systematically. We fill this gap by proposing a formalized technique that takes into consideration data reuse, limited lifetime of the arrays of an application and application specific prefetching opportunities, and performs a thorough tradeoff exploration for different memory layer sizes. This technique has been implemented on a prototype tool, which was tested successfully using nine real-life applications of industrial relevance. Following this approach we have able to reduce execution time up to 60%, and energy consumption up to 70%.

Description of our work

We overcome the performance and energy bottleneck that memory imposes at the system by developing the Memory Allocation and Layer Assignment [START_REF] Brockmeyer | Layer assignment techniques for low energy in multi-layered memory organizations[END_REF] with time extensions (MHLA with TE) technique. Our technique targets at the data reuse opportunities, found in most multimedia and image processing applications, array in-place optimizations and optimum memory transfer scheduling. By exploiting data reuse, a part of an array is copied from one layer to a lower layer, closer to the processor. As a result, energy and performance can be improved since most accesses take place on the smaller copy and not on the large more energy consuming and slower higher memory layer. TE enables us to increase the performance of a system, with minimal cost. They enable the selective prefetching of copy can-didates from off-chip memory layers to on-chip memory layers, exploiting the lifetime information of copies. Time extensions, need the support of a memory transfer engine (like DMA engine or data mover) that allows simultaneous the CPU to continue processing data and the engine to copy off-chip data to on-chip layers. The goal of time extensions is to hide as much as possible the cycles required in accessing off-chip memory, respecting data dependencies and on-chip size requirements. In case that our architecture does not support a memory transfer engine, TE are not applicable.

MHLA and Time Extensions

The MHLA technique with time extensions has been implemented in a tool with the same name (MHLA), which is able to find all the optimal trade-off points, given some architecture specific constraints and models. The exploration flow can be divided into two distinct steps: A selection and assignment step and a time extension step. Every step has a number of substeps. When we first presented the MHLA technique [START_REF] Brockmeyer | Layer assignment techniques for low energy in multi-layered memory organizations[END_REF], it did not consider TE. We present a new version of our MHLA technique with the major contribution the support of TE. In our technique TE is the second step in MHLA.

Time extensions are done in an iterative process (Figure 1). We examine every DMA Block Transfer (BT) and we try to schedule earlier the initiating of the DMA, obeying dependencies and on-chip memory size.

We iterate over the list of BTs in the greedy order and try to perform prefetching. We evaluate whether the increase in lifetime of the copy candidate, due to the extension of the DMA one loop before, is valid or not. If the increase in lifetime increases on-chip memory size over the user-defined on-chip memory constraint, then this extension is not valid and no further actions are performed for this BT. Otherwise, TE for one loop before is valid for this BT.

Experimental Results

The usefulness and feasibility of the MHLA technique is demonstrated using nine real life applications of motion estimation, video encoding, image and audio processing domain. Although, we only consider single threaded applications, we plan to extend our technique to multiple tasks with multiple threads.

MHLA technique performs exploration in two steps. After deciding and placing on memory layers, arrays and copies the step of time extensions is applied. On the applications considered, the first step boost performance from 40% to 60% compared to the out of the box code (Figure 2) for specific memory sizes. An optimum memory allocation and assignment can also reduce energy consumption significantly up to 70% (Figure 3).

The MHLA second step tries to perform prefetching on the copies that have been decided in the previous step. This step can boost performance of up 33%, if there are a lot of processing loops that can hide prefetching block transfers. This step pushes performance towards the ideal case where ever block transfer can be hidden from the processor (0 wait cycles block transfer time) (Figure 2). Energy consumption in both steps remains the same, because in our models we only consider accesses to the memory hierar- chy. Though, it is expected that the reduced execution time will result in lower energy consumption.

Figure 1 .

 1 Figure 1. The TE step performs an application specific prefetching.

Figure 2 .

 2 Figure 2. MHLA improves performance up to 60%. MHLA with TE can boost performance even more.

Figure 3 .

 3 Figure 3. In addition to performance improvements, MHLA technique benefits energy consumption as well.

* This work was sponsored by a scholarship from Public Benefit Foundation of Alexander S. Onassis and from Marie Curie Host Fellowship project HPMT-CT-2000-00031.

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Conclusions

We have addressed the performance/energy consumption bottleneck problem by developing a a formalized technique of memory allocation and assignment. This technique, exploits data-reuse opportunities, in-place optimizations and application specific prefetching, which should be used during the early system design steps. This technique has been implemented in a prototype tool, which is part of the ATOMIUM framework, and has allowed us to do fast, accurate and automatic exploration of nine real-life applications. We found significant performance and energy consumption gains on every application, illustrating the importance of our technique in the early design phases.


Peng Li

Frank Liu

Xin Li

Lawrence T Pileggi

Sani R Nassif

Modeling Interconnect Variability Using Efficient Parametric Model Order Reduction

Assessing IC manufacturing process fluctuations and their impacts on IC interconnect performance has become unavoidable for modern DSM designs. However, the construction of parametric interconnect models is often hampered by the rapid increase in computational cost and model complexity. In this paper we present an efficient yet accurate parametric model order reduction algorithm for addressing the variability of IC interconnect performance. The efficiency of the approach lies in a novel combination of low-rank matrix approximation and multi-parameter moment matching. The complexity of the proposed parametric model order reduction is as low as that of a standard Krylov subspace method when applied to a nominal system. Under the projection-based framework, our algorithm also preserves the passivity of the resulting parametric models.

Introduction

The every-increasing variability of modern IC manufacturing process has introduced significant variations in circuit performances. Accordingly, assessment of these performance variations has become an integral part of analysis and optimization of modern VLSI designs. In the past decade, efficiency in IC interconnects analysis has been facilitated by powerful model order reduction techniques [START_REF] Pillage | Asymptotic waveform evaluation for timing analysis[END_REF]- [START_REF] Odabasioglu | PRIMA: passive reducedorder interconnect macromodeling algorithm[END_REF]. Two classes of model order reduction techniques have gained popularity: moment-matching based method (e.g. AWE [START_REF] Pillage | Asymptotic waveform evaluation for timing analysis[END_REF] and Krylov subspace methods [START_REF] Feldmann | Efficient linear circuit analysis by Padé approximation via the Lanczos process[END_REF]- [START_REF] Odabasioglu | PRIMA: passive reducedorder interconnect macromodeling algorithm[END_REF]) and control-theoretic approaches (e.g. truncated balanced realization (TBR) [START_REF] Moore | Principle component analysis in linear systems: controllability, observability, and model reduction[END_REF][8] [START_REF] Wong | Passivity-preserving model reduction via a computationally efficient project-and-balance scheme[END_REF]). The first class of methods are very attractive in terms of computational cost while the methods that fall into the second class tend to be more accurate, but suffer from a dramatic increase in computational cost. The high cost associated with the latter often preclude these methods from being directly applied to large practical problems.

Variational modeling of interconnects can be classified similarly: moment-matching (Krylov subspace) based [START_REF] Liu | Model order reduction of RC(L) interconnect including variational analysis[END_REF] vs. TBR-based [START_REF] Heydari | Model reduction of variable-geometry interconnects using variational spectrally-weighted balanced truncation[END_REF][12] [START_REF] Wang | A linear fractional transform (LFT) based model for interconnect parametric uncertainty[END_REF]. In [START_REF] Liu | Model order reduction of RC(L) interconnect including variational analysis[END_REF], Liu et. al. proposed to construct variational models by fitting the projection vectors over the samples taken in the variational parameter space. In [START_REF] Heydari | Model reduction of variable-geometry interconnects using variational spectrally-weighted balanced truncation[END_REF], Heydari et. al. performed TBR analysis on perturbed systems. Additionally, a linear fractional transform (LFT) formulation based TBR-like approach was proposed in [START_REF] Beck | Model reduction of multidimensional and uncertain systems[END_REF] [START_REF] Wang | A linear fractional transform (LFT) based model for interconnect parametric uncertainty[END_REF]. Almost the same comparison can be drawn for these variational methods: moment matching based approaches are computationally more efficient while TBR methods tend to be more accurate but with a significantly higher cost.

Nevertheless, in comparison to modeling techniques employed for the nominal system, the inclusion of variational effects into the model order reduction implies a significant increase in computational cost and/or model complexity. This can be true even for the more efficient moment matching methods. Hence, in this paper we focus on developing more efficient Krylov subspace methods.

In this paper, we first show that in principle, variational modeling of interconnects can be cast as a general multiparameter Krylov-based moment matching problem [START_REF] Daniel | A multiparameter moment-matching model-reduction approach for generating geometrically parameterized interconnect performance models[END_REF]. We then point out that for the particular problem of modeling process variations this is not an entirely suitable choice due to its limitation in model complexity. As an improvement, we show that exploiting multi-point expansion under the same Krylov framework can lead to more compact interconnect models but with an increase in computational cost. Furthermore, we propose a novel combination of lowrank matrix approximation and multi-parameter moment matching scheme and show that we can maintain the same low cost of standard Krylov methods while significantly improve model complexity. The new approach allows us to achieve a higher moment-matching order for multiparameter moments than what was possible before, thereby significantly improving the model accuracy. One further advantage of the approach is that the passivity of reduced parametric models can be easily guaranteed. We demonstrate the proposed technique on several variational interconnect modeling examples.

Prior Work

An interconnect network can be described by the following MNA (modified nodal analysis) formulation:

x L y Bu Gx x C T = + - = , D (1)
where G , C , B and L are system matrices, m R u ∈ and m R y ∈ denote port voltages and currents, n R x∈ is the vector consisting of nodal voltages and branch currents for voltage sources and inductors. The PRIMA algorithm of [START_REF] Odabasioglu | PRIMA: passive reducedorder interconnect macromodeling algorithm[END_REF] computes an orthonormal basis

q n R V × ∈ for a Krylov subspace of C G A 1 - - ≡
, which is spanned by several block moments of the system transfer function. The reduced order model of ( 1) is obtained by computing the following matrices of smaller dimensions:

L V L B V B CV V C GV V G T T T T = = = = , , , ~ .
(2) To model the variational effects of interconnects, [START_REF] Liu | Model order reduction of RC(L) interconnect including variational analysis[END_REF] uses the first order expansions to express the system conductance and capacitance matrices:

2 2 1 1 0 2 1 2 2 1 1 0 2 1 ) , ( , ) , ( C p C p C p p C G p G p G p p G + + = + + = (3)
G , 1 C , 2 G , 2
C the sensitivity matrices w.r.t variational parameters. When the circuit parameters vary due to process fluctuations, the projection matrix used in PRIMA also undergoes small variations. The projection matrix is therefore expanded using a Taylor series [START_REF] Liu | Model order reduction of RC(L) interconnect including variational analysis[END_REF]:

2 2 22 2 12 2 1 21 1 11 0 2 1 ) , ( p V p V p V p V V p p V + + + + =
(4) To determine the coefficient matrices in the above equation, samples are taken in the variational parameter space and the PRIMA algorithm is applied to the resulting perturbed systems. Based on these sampled projection matrices, coefficient matrices in (4) are determined by solving a set of linear equations. Inserting (4) into (2), a variational reduced order model parametrizable in s p i ' is obtained.

Multi-Parameter Moment Matching

In this section, we first point out that parametric reduced order modeling can be more generally cast as a multiparameter Krylov projection based moment matching problem. In what follows, we show the benefit of exploiting multi-point expansions under this context.

Single-point expansion

A Krylov projection based multi-parameter moment matching approach was adopted in [START_REF] Daniel | A multiparameter moment-matching model-reduction approach for generating geometrically parameterized interconnect performance models[END_REF] to generate parameterized interconnect performance models. We shall refer to the approach in [START_REF] Daniel | A multiparameter moment-matching model-reduction approach for generating geometrically parameterized interconnect performance models[END_REF] as a single-point expansion method since the system transfer function is expanded at only one point in the joint space of s and s p i ' .

+ + + + = + + + + = 2 2 1 1 0 2 1 2 2 1 1 0 2 1 ) , , , ( ) , , , ( (5)
We use the notation

p p p n n n L B C C G G C G } , , , , , , , , , { 1 1
0 0 to specify a p n -parameter dynamic system as in (5). The transfer function of the system can be written as

[ ] B G C G sp G G p C sG I B p p p sC p p p G p p p s X p p p p p n i n i i i i i n n n 1 0 1 1 1 1 0 1 0 0 1 0 1 2 1 2 1 2 1 ) , , , ( ) , , , ( ) , , , , ( -
- = = - - - -       + + + = + = ∑ ∑ (6)
k n k k k k k k k k k k k k k k k k k k k k k k k n i n i i i i i n p p s M B G C G sp G G p C sG p p p s X ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∞ = = - = - - = - - - = ∞ = - = = - - - - =       - - - = 0 0 0 0 0 1 , , , , 0 1 0 1 1 1 0 1 0 0 1 0 2 1 1 1 2 1 1 1 1 ) , , , , ( (7)
In ( 7)

)

( 1 , , , , 1 k k k k M p p n s n s k k k k = + + is a k-th order multi- parameter moment corresponding to the coefficient of the term p n p s k n k k k p p p s 2 1 2 1
. Consider an orthonormal basis

) ( n q R V q n < ∈ ×
of the subspace spanned by the multiparameter moments with an order equal to or less than k

} , , , , { ) ( , , , 0 , 1 , 1 0 , 0 , 0 1 p n s k k k k M M M span V colspan = (8)
We can employ V as a projection matrix to compute a reduced model; i.e., each system matrix of the full model including sensitivity matrices is reduced using matrix transformations of (2). For instance, the sensitivity of

k p G is reduced to V G V G k k p T p = ~.
It can be shown that the resulting reduced model matches up to the k-th order multi-parameter moments of the original transfer function.

Inefficiency of single-point expansion

Although all the parameters are treated identically in the above multi-parameter moment matching scheme, there does exist an asymmetry between these parameters that arise practically in our manufacturing-variation-aware interconnect modeling problem. While variational parameters fluctuate in a small range around the nominal values due to process variations, the "variation" in the frequency variable s has to be considered over a wider range in order to capture the complete system frequency response.

It is not surprising that the reduced order model size of the above approach grows rapidly with the number of parameter and the orders of moment matching. With more careful analysis, it is quite interesting to see that the model size actually depends not so much on the moment matching order of individual parameters but the highest moment matching order with respect to any parameter. This is because a high moment matching order of any parameter will introduce a large number of cross terms in the multi-parameter expansion, thus lead to a large reduced model size. In other words, for interconnect variational modeling, choosing a low moment matching order for variational parameters (due to the small range of variation) does not significantly reduce the model size if the order for variable s remains high. As a result, the single-point expansion produces an expensive reduced model even when the circuit is only perturbed within a small range about the nominal condition.

Multi-point expansion

Let us consider the multi-point expansion in Fig. 1, where multiple samples are taken within the specified ranges in the variational parameter space. Let us assume that s n samples are taken and for each perturbed system

T n i n i i i p p p p p p p p P ] , , , [ , 2 , 2 1 , 1 = = = =
, a projection matrix i V is computed to match the k moments of s using any standard Krylov method. Finally, a parametric model can be produced by forming an orthonormal basis V of the

combined projection vectors } , , , { 2 1 s n V V V colspan V
= and employing V as the final projection matrix. The resulting model preserves the first k moments of s at the above s n sample points. In some sense, the reduced model approximates the full model at the sample points using a standard model order reduction technique, and then interpolates implicitly between these samples in the variational parameter space. In addition to other benefits of a multi-point method, it can be used to reduce the model size. To see this, let us assume that we want to match multi-parameter transfer function moments corresponding to the following terms:

i k i i i i k i i i p s p s sp p p s p s p s p s , , , , , , , , , 2 0 0 2 0 1 0 0
, where the expansion order for some variational parameter i p is the 1st order. To do so, using a single-point method will lead to a reduced order model of size

m k k ) 1 ( 2 + +
, where m is the number of ports. In a quite equivalent way, one can employ a multi-point method to match the first . A similar situation is encountered under the context of nonlinear system model order reduction [START_REF] Li | NORM: Compact model order reduction of weakly nonlinear systems[END_REF].

The aforementioned multi-point method bears resemblance of the approach of [START_REF] Liu | Model order reduction of RC(L) interconnect including variational analysis[END_REF]. The major difference is that in [START_REF] Liu | Model order reduction of RC(L) interconnect including variational analysis[END_REF] interpolation between the sample points is done by direct fitting while in this approach it is accomplished implicitly via projection. Sometimes it is observed that the projection matrix is sensitive w.r.t variational parameters thus making a direct fitting less robust. Under these cases, multi-point expansion might be a more robust choice.

Low-Rank Approximation Based Multi-Parameter Moment Matching 4.1 The proposed approach

Although sampling in variational parameter space can lead to an improvement in model complexity, the computational cost and/or model size can still be high when several variational parameters are considered simultaneously. For instance, if we sample a four-dimensional parameter space by taking three samples per axis, a total of 81 sample points will result that requires the same number of matrix factorizations. To see the interactions between different parameters, we rewrite the transfer function expansion in (7)

( ) ( ) 0 1 1 1 1 , , 1 0 1 0 0 2 1 , ) , , , , ( R C G G G A p p s I p p p s X k n i n j k l c g a c j g i a c j g i c a n p p l l l l l l l l l l p                 - - + = ∑∑∑ ∑ ∞ = = = = - - + (9) where 0 1 0 0 C G A - - = , B G R 1 0 0 - = , } 1 , 0 { , , ∈ l l l c g a and 1 = + + l l l c g a . We refer to each i G G 1 0 - , i C G 1 0 - , p n i ≤ ≤ 1 as
a generalized sensitivity matrix with respect to 0 G .

Consider one possible term in the summation of ( 9)

( ) 0 0 1 0 0 R A G G A p s F q i t i q t - + - = (10)
To match the moment expression in (10) more efficiently, we first seek the optimal 2-norm rank-

i G k ) ( n k i G << approximation of the generalized sensitivity matrix i G G 1 0 - using SVD (singular value decomposition) [16] ] , , [ ], , , [ ˆ, ˆ1 1 1 1 1 0 i G i i G i G i i i i G k G k k G T G G k j T j j j i v v V u u U V U v u G G = = ⋅ = ≈ - ∑ = - σ σ σ ( 11
)
where

i G j j j k j v u ≤ ≤ 1 , , , σ are the i G k largest singular values
and the corresponding left and right singular value vectors, respectively. Substituting ( 11) into (10) yields 0 0

1 1 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 ) )( ( ˆ) ( ) )( ( ˆR C G C G V U A p s R C G C G C G V U A p s R A V U A p s F q T G G t i q t T G G t i q t q T G G t i q t i i i i i i - - - + - - - + + - - = - - - = ≅ (12) Now defining i i G T G V G V 0- - = and T T T C G A 0 0 0 - - = , we have [ ] [ ] 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 1 0 0 1 0 0 ) ) ( ˆ) ( ) ( ˆ) )( ( ˆR C V A U A p s R C V G C G U A p s R C G C G V U A p s F T G q T G t i q t T G T q T T G t i q t q T G G t i q t i i i i i i - + - - - + - - - + = - - = - - ≅ (13)
There are two important implications due to the approximation in [START_REF] Wang | A linear fractional transform (LFT) based model for interconnect parametric uncertainty[END_REF]. First, it suggests that there are two useful Krylov subspaces, one w.r.t 0 A and the other w.r.t T A 0 that can be used for matching the corresponding moment:

} , , { ) , , ( } , , , { ) 1 , , ( 1 0 0 0 0 0 0 i i i i i i i i G q T G T G G T G t G G G V A V A V colspan q V A Kr U A U A U colspan t U A Kr - = = + (14)
Secondly, the use of low-rank matrix approximation for generalized sensitivity matrices has decoupled the Krylov subspace construction w.r.t various parameters, i.e. these subspaces can be computed independently from each other and combined at a later time to build the overall projection matrix. Although we have only shown this decoupling between frequency variable s and parameter i p in the simple case of (10) due to the space limitation, the same property remains even for an arbitrary combination of multiple parameters. This comes from the fact that we include the Krylov subspaces associated with all the lowrank approximation vectors in the projection. The immediate q terms benefit of this result is that the size of the reduce order model can be significantly reduced as the cross terms of multi-parameter moments are no longer an issue.

It is worth mentioning that a similar reduction in model size can be achieved if low-rank approximations are applied to sensitivity matrices instead of generalized sensitivity matrices. However, this choice will incur a larger error in approximating the true transfer function. We have observed that approximating the generalized sensitivity matrices work much better in practice due to their stronger connection to moments, as evident in [START_REF] Li | NORM: Compact model order reduction of weakly nonlinear systems[END_REF]. We outline the proposed lowrank approximation based single-point multi-parameter moment matching scheme in algorithm 1 (Fig. 2). We summarize the moment matching property of algorithm 1 in theorem 1. Due to the space limitation, we omit the proof of the theorem. Theorem 1 (moment matching property of algorithm 1):

For the approximate parametric system p p p n n n L B C C G G C G } , , , , , , , , , { 1 1 0 0 m m , where T G G i i i V U G G 0 - = , T C C i i i V U G C 0 - = , p n i ≤ ≤ 1
, the reduced order model using the projection matrix computed in algorithm 1 preserves the transfer function moments with respect to all parameters up to order k.

As in Theorem 1, the projection V computed in Algorithm 1 corresponds to a nearby system that is based on the lowrank approximation of generalized sensitivity matrices. In other words, we compute an approximate projection matrix for the original parametric system and employ it for parametric model order reduction. It should be noted that the same moment matching property of theorem 1 can be achieved by neglecting the Krylov subspaces with respect to V ˆto the projection matrix. This simplification can reduce the model size approximately by a factor of two. Nevertheless, once the projection matrix is computed we will apply matrix reduction on the sensitivity matrices but not their low-rank approximations as in step 4 of the algorithm. Therefore, incorporating the useful Krylov subspaces of T A 0 improves the accuracy of the model order reduction. The congruence transforms employed in step 4 of the algorithm also implies that the passivity of the reduced model will be guaranteed if the original parametric model is passive.

Computational cost and model complexity

One appealing property of Algorithm 1 is that the dominant cost of the parametric model order reduction is the same as that of a Krylov subspace method when applied to a deterministic system. The dominant cost is again one-time factorization of a sparse matrix, 0 G .

Although the generalized sensitivity matrices are dense, the computation of a few dominant singular values/vectors for them (step 1 of algorithm 1) can be efficiently facilitated by the matrix-implicit nature of iterative sparse SVD techniques [START_REF] Berry | Large scale singular value computations[END_REF] [START_REF] Larsen | Lanczos bidiagonalization with partial reorthogonalization[END_REF]. For instance, a low-rank approximation of

T i G G 1 0 - -
can be efficiently done using a few subspace iterations wherein the dense generalized sensitivity matrix is not explicitly required but only its matrix-vector products. The latter can be provided efficiently by reusing the same factorization of 0 G . Furthermore, the computation of Krylov subspaces with respect to T A 0 can be achieved by the same sharing. Notice that if the LU factorization of 0 G is

g g U L G = 0 , then T g T g T L U G = 0 . The matrix-vector product x A y T ⋅ = 0 can be achieved by solving x C y G T T 0 0 = - or x C y L U T T g T g 0 = -
. Therefore, the cost of the algorithm is linear in both the moment matching order k and the number of variational parameters p n , and almost linear in the number of circuit nodes.

For a k-th order multi-parameter moment matching, the size of reduced order model is

k m n k O p svd ) 4 ( + , where svd k
is the rank of the SVD approximation, and m the number of ports (using the simplification pointed out in section 4.1

Algorithm 1:

1. Compute a low-rank approximation for each generalized relative sensitivity matrix:

T G G i i i V U G G 1 0 ≈ -- , T C C i i i V U C G 1 0 ≈ -- , p n i ≤ ≤ 1

Compute an orthonormal basis of the following

Krylov subspaces:

2.1. } , , , { 0 0 0 0 0 0 R A R A R colspan V k m = where B G R C G A 1 0 0 0 1 0 0 , - - = - = 2.2. For each p n i i ≤ ≤ 1 , , i i G T G V G V 0- - = , i i C T C V G V 0- - = } , , , { 0 0 1 , i i i i G k G G G U A U A U colspan V m = } , , , { 1 0 0 2 , i i i i G k T G T G G V A V A V colspan V - = m } , , , { 1 0 0 1 , i i i i C k C C C U A U A U colspan V - = m } , , , { 2 0 0 2 , i i i i C k T C T C C V A V A V colspan V - = m where T T T C G A 0 0 0 - - =

Compute an othonormal basis V for the combined

subspaces: leads to a model size of

} , , , , , , , , , { 2 , 1 , 2 , 1 , 2 , 1 , 2 , 1 , 0 1
1 1 1 p n p n p n p n C C G G C C G G V V V V V V V V V colspan V m = 4. Construct the final reduced order model: L V L B V B V C V C V G V G T T T T = = = = , , , ~0 0 0 0 , V G V G i T i = ~, V C V C i T i = ~, for p n i ≤ ≤ 1
k m n k O p svd ) 2 ( + ).
In practice, we have observed that a rank-one approximation is usually sufficient to provide a good accuracy. This model complexity is in contrast to that of a multi-point expansion approach,

) ( mk c O p n
, in which c samples per parameter are taken in the p n -dimensional variational parameter space.

Numerical Results

An RC network

We consider an RC network of 767 circuit unknowns subjected to two independent variational sources. We randomly vary the RC values of the circuit, and then extract the sensitivity matrices w.r.t. these two variational sources correspondingly. Finally, we applied the proposed algorithm to construct a variational reduced order model of size 37 (number of states). As outlined in Fig. 2, this model (approximately) matches up to 4 th order multi-parameter moments w.r.t all parameters of the full model. For comparison, by taking 8 samples in the variational parameter space, we compute a 40-state multi-point expansion based model that matches up to 4 th order moments w.r.t s at each sample point. We evaluate the model accuracy for a perturbed network obtained by injecting up to 70% parametric variations into the nominal system. In Fig. 3, the transfer function from the voltage input to an observation node are plotted for five models: the nominal system, the perturbed system and three reduced models of the perturbed system. In addition to the aforementioned two reduced order models, we also include the third reduced order model in Fig. 3., which is obtained by projecting the variational model using a projection matrix obtained by applying PRIMA algorithm to the nominal system (matching 8 moments of s). It can be clearly seen from Fig. 3 that while the model based on the nominal projection matrix fails to capture the variation in transfer function, the other two reduced models become almost not indistinguishable with the perturbed full model in the plot.

A 4-port RLC network

Next, a two-bit bus is modeled as a coupled 4-port RLC network, where each line consists of 180 RLC segments. The size of MNA formulation for the bus is 1086. We follow the same procedure as in the previous example to generate parametric variations for two independent variational sources. Again, we consider three variational reduced order models. The first one has a size of 52, is obtained by employing the nominal projection matrix. The second model has a size of 144 and is computed by applying the proposed approach to (approximately) match the moments of all parameters (including cross-terms) up to 12 th order. Among the moments matched, 52 are multi-port moments w.r.t s. The last reduced model is generated by taking three samples in the variational parameter space and applying the multi-point expansion. This model has a size of 156 and matches 52 multi-port moments w.r.t. s at each sample point. Various models are examined by evaluating a perturbed system subjected to a maximum of 30% parametric variation in Fig. 4. In the figure, we plot port admittance 11 Y based on the nominal full model, perturbed full model and three reduced order models. As shown in the plot, , the frequency domain response of a RLC circuit becomes more sensitive to parametric variations. It is also evident that for this case that building a variational model based upon only the information of the nominal system is far from adequate. However, by incorporating the variational information effectively into the model order reduction, the low-rank approximation based model is able to capture the variations in circuit response fairly accurately. In comparison, the larger multi-point-method model is not as accurate and the domain computational cost is three times larger.

Clock Trees

For more realistic examples, we consider two industrial RC networks, RCNetA and RCNetB. They are portions of a clock tree, and routed on three metal layers: M5, M6 and M7. RCNetA has 78 nodes while RCNetB 333. We consider three independent metal line width variations on these metal layers. In this case, the sensitivity matrices w.r.t metal line width variations are obtained by performing multiple parasitic extractions.

For RCNetA, we applied the low-rank based parametric model reduction to compute a reduced order model of size 29 while matching the moments of s to the 4 th order and the rest of multi-parameter moments to the 2 nd order. We independently vary the three metal line widths up to 30% ( σ 3 variations) of the nominal values according to the normal distribution. For each perturbed circuit instance, we compare the reduced order model with the full model in terms of the 5 most dominant poles. The error distribution in these poles across all the instances is plotted in Fig. 5 (left plot). For this relatively small example, the error in these dominant poles is completely negligible. In Fig. 5 (right plot), we show the error in the most dominant pole as a function of M5 and M6 metal line widths (within -30% to 30% of their nominal values). This again confirms the accuracy of the reduced order model. The same steps are taken to examine the model accuracy for RCNetB. We first compute a reduced order parametric model of size 40 while matching all the multi-parameter moments to the 3 rd order. Then, we independently vary the three metal line widths up to 30% ( σ 3 variations) of the nominal values using the normal distribution and plot the error in the 5 most dominant poles across all the instances in Fig. 6 (left plot). The maximum error out of 1000 poles is less than 0.12%. Next, we examine how the error in the most dominant pole changes as we vary the metal line widths of M5 and M6 in between -30% and 30%. As shown in Fig. 6 (right plot), the reduced order model is very accurate as the largest error is less than 0.3%.

Conclusion

We have shown that multi-point expansion in variational parameter space can be exploited to produce more compact parametric interconnect models than its single-point counterpart. Furthermore, we show that by combining lowrank matrix approximation and multi-parameter moment matching even more efficient parametric models can be constructed. Additionally, these variation-aware models are obtained at an almost linear cost and preserve the passivity. Relative Error of The Dominant Pole (%) M6 Variation (%)

-30% M5 Variation -15% M5 Variation 0% M5 Variation 15% M5 Variation 30% M5 Variation

Cp and 2 p

 2 are conductance and capacitance matrices under the nominal conditions, 1 p and 2 p are parameters used to model variational sources such as the metal line width or thickness variations. For convenience, we refer to 1 as variational parameters, and 1

Fig 1 .

 1 Fig 1. Multi-point expansion in the process parameter space.

Fig 2 .

 2 Fig 2. Low-rank approximation based single-point multi-parameter moment matching.

  Fig 3. A RC net for which the RC elements are subjected to a maximum 80% parametric variation.

Fig 4 . A 4 -

 44 Fig 4. A 4-port coupled RLC network. The perturbed system is subjected to a maximum 30% parametric variation.

Fig 5 .

 5 Fig 5. Model accuracy for RCNetA. Left: error distribution in the 5 most dominant poles; Right: the error in the most dominant pole as a function of metal width variations.

Fig 6 .

 6 Fig 6. Model accuracy for RCNetB. Left: error distribution in the 5 most dominant poles; Right: the error in the most dominant pole as a function of metal width variations. 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 50 100 150 200 250 300 350 400

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE