R Jha

D Chablat
email: damien.chablat@irccyn.ec-nantes.fr

F Rouillier
email: fabrice.rouillier@inria.fr

G Moroz

Workspace and Singularity analysis of a Delta like family robot

Keywords: Delta-like robot, Cylindrical algebraic decomposition, Workspace, Gröbner basis, Parallel robot

Workspace and joint space analysis are essential steps in describing the task and designing the control loop of the robot, respectively. This paper presents the descriptive analysis of a family of delta-like parallel robots by using algebraic tools to induce an estimation about the complexity in representing the singularities in the workspace and the joint space. A Gröbner based elimination is used to compute the singularities of the manipulator and a Cylindrical Algebraic Decomposition algorithm is used to study the workspace and the joint space. From these algebraic objects, we propose some certified three dimensional plotting describing the the shape of workspace and of the joint space which will help the engineers or researchers to decide the most suited configuration of the manipulator they should use for a given task. Also, the different parameters associated with the complexity of the serial and parallel singularities are tabulated, which further enhance the selection of the different configuration of the manipulator by comparing the complexity of the singularity equations.

Introduction

The workspace can be defined as the volume of space or the complete set of poses which the end-effector of the manipulator can reach. Many researchers published several works on the problem of computing these complete sets for robot kinematics. Based on the early studies [START_REF] Soni | An algorithm for the workspace of a general n-R robot[END_REF][START_REF] Hansen | Generation and evaluation of the workspace of a manipulator[END_REF], several methods for workspace determination have been proposed, but many of them are applicable only for a particular class of robots. The workspace of parallel robots mainly depends on the actuated joint variables, the range of motion of the joints and the mechanical interferences between the bodies of mechanism. There are different techniques based on geometric [START_REF] Gosselin | Determination of the workspace of 6-DOF parallel manipulators[END_REF][START_REF] Merlet | Geometrical determination of the workspace of a constrained parallel manipulator[END_REF], discretization [START_REF] Castelli | A fairly general algorithm to evaluate workspace characteristics of serial and parallel manipulators[END_REF][START_REF] Bonev | A new approach to orientation workspace analysis of 6-DOF parallel manipulators[END_REF][START_REF] Chablat | An Interval Analysis Based Study for the Design and the Comparison of 3-DOF Parallel Kinematic Machines[END_REF], and algebraic methods [START_REF] Zein | An exhaustive study of the workspace topologies of all 3R orthogonal manipulators with geometric simplifications[END_REF][START_REF] Ottaviano | Identification of the workspace boundary of a general 3-R manipulator[END_REF][START_REF] Chablat | Non-singular assembly mode changing trajectories in the workspace for the 3-RPS parallel robot[END_REF][START_REF] Chablat | Uniqueness domains and non singular assembly mode changing trajectories[END_REF] which can be used to compute the workspace of parallel robot. The main advantage of the geometric approach is that, it establish the nature of the boundary of the workspace [START_REF] Siciliano | Springer handbook of robotics[END_REF]. Also it allows to compute the surface and volume of the workspace while being very efficient in terms of storage space, but when the rotational motion is included, it becomes less efficient. Interval analysis based methods can be used to compute the workspace but the computation time depends on the complexity of the robot and the requested accuracy [START_REF] Chablat | An Interval Analysis Based Study for the Design and the Comparison of 3-DOF Parallel Kinematic Machines[END_REF]. Discretization methods are usually less complicated and can easily take into account all kinematic constraints, but they require more space and computation time for higher resolutions. The majority of numerical methods which is used to determine the workspace of parallel manipulators includes the discretization of the pose parameters for computing workspace boundaries [START_REF] Bonev | A new approach to orientation workspace analysis of 6-DOF parallel manipulators[END_REF]. There are other approaches, which are based on optimization algorithms [START_REF] Snyman | An optimization approach to the determination of the boundaries of manipulator workspaces[END_REF] for fully serial or parallel manipulators, analytic methods for symmetrical spherical mechanisms [START_REF] Bonev | Analytical determination of the workspace of symmetrical spherical parallel mechanisms[END_REF]. In [START_REF] Bohigas | A complete method for workspace boundary determination on general structure manipulators[END_REF] a method for computing the workspace boundary for manipulators with a general structure is proposed, which uses a branch-and-prune technique to isolate a set of output singularities, and then classifies the points on such set according to whether they correspond to motion impediments in the workspace. A cylindrical algebraic decomposition (CAD) based method is illustrated in [START_REF] Chablat | Non-singular assembly mode changing trajectories in the workspace for the 3-RPS parallel robot[END_REF][START_REF] Chablat | Workspace and joint space analysis of the 3-RPS parallel robot[END_REF], which is used to model the workspace and joint space for the 3 RPS parallel robot. This paper presents the results which are obtained by applying algebraic methods for the workspace and joint space analysis of a family of delta-like robot including complexity information for representing the singularities in the workspace and the joint space. The CAD algorithm is used to study the workspace and joint space, and a Gröbner based elimination process is used to compute the parallel and serial singularities of the manipulator. The structure of the paper is as follows. Section 2 describes the architecture of the manipulator, including kinematic equation and joint constraints associated with the manipulators. Section 3 discusses the computation of parallel as well as serial singularities and their projections in workspace and jointspace. Section 4 presents a comparative study on the shape of the workspace of different delta-like robots. Section 5 finally concludes the paper.

Manipulators Architecture

The manipulator architecture under the study is a three degree of freedom parallel mechanism which consists of three identical legs, the different arrangements of these legs give rise to family of delta like robot. Several types of delta-like robot were studied, few of them are Orthoglide [START_REF] Chablat | An Interval Analysis Based Study for the Design and the Comparison of 3-DOF Parallel Kinematic Machines[END_REF][START_REF] Pashkevich | Kinematics and workspace analysis of a three-axis parallel manipulator: the Orthoglide[END_REF], Hybridglide, Triaglide [START_REF] Hebsacker | Hexaglide 6 DOF and Triaglide 3 DOF Parallel Manipulators[END_REF] and UraneSX [START_REF] Chablat | An Interval Analysis Based Study for the Design and the Comparison of 3-DOF Parallel Kinematic Machines[END_REF]. The position vectors of end points of i th leg are P j and B j i , also A j i and B j i are the position vectors of end points of i th actuator where j represents the manipulator type from the family of delta like robot ( j = 1 -Orthoglide, 2 -Hybridglide, 3 -Triaglide, 4 -UraneSX). ρ ρ ρ ρ ρ ρ j i represents the prismatic joint variables whereas P j represents the position vector of the tool center point which is shown in Eq. ( 1). The kinematic equations of the family of delta like robot can be generalized as

||A j i B j i || = ρ j i P = [x y z] T with i = 1, 2, 3 and j = 1, 2, 3, 4 (1)
||P -B j i || = L j i .
All the computations and analysis are done for L j i = L and by imposing the following constraints on joint variables. Without joint limits, all the robots admit two assembly modes and eight working modes.

0 < ρ 1 < 2L 0 < ρ 2 < 2L 0 < ρ 3 < 2L (2)
The Orthoglide mechanism is driven by three actuated orthogonal prismatic joints.

A simpler virtual model can be defined for the Orthoglide, which consists of three bar links connected by the revolute joints to the tool center point on one side and to the corresponding prismatic joint at another side. Several assembly modes of these robots depends upon the solutions of direct kinematic problem is shown in Fig. 1(a). The point P i represents the pose of corresponding robot. However more than one value of i for the point P i shows the multiple solutions for the DKP. The constraint equations for the Orthoglide are

(x -ρ 1 ) 2 +y 2 +z 2 =L 2 ; x 2 +(y -ρ 2 ) 2 +z 2 =L 2 ; x 2 +y 2 +(z -ρ 3 ) 2 =L 2
The Hybridglide mechanism consists of three actuated prismatic joints, in which two actuators are placed parallel and third one perpendicular to others two. Also the three bar links connected by spherical joints to the tool center point on one side and to the corresponding prismatic joint at another side. Several assembly modes of these robots depends upon the solutions of direct kinematic problem is shown in Fig. 1(b) and following is the constraint equations:

(x -1) 2 +(y-ρ 1 ) 2 +z 2 = L 2 ; (x +1) 2 +(y-ρ 2 ) 2 +z 2 = L 2 ; x 2 +y 2 +(z-ρ 3 ) 2 = L 2
The Triaglide manipulator is driven by three actuated prismatic joints, in which all the three actuators parallel to each other and placed in the same plane. The architecture of the Triaglide is shown in Fig. 1(c) and the constraint equations are defined as:

(x -1) 2 + (y -ρ 1 ) 2 + z 2 =L 2 ; (x + 1) 2 + (y -ρ 2 ) 2 + z 2 =L 2 ; x 2 + (y -ρ 3 ) 2 + z 2 =L 2
The UraneSX is similar to triaglide, but instead of three actuators in the same plane, it is placed in different planes. The architecture of the Triaglide is shown in Fig. 1(d) and the constraint equations are as follows:

(x -1) 2 + y 2 + (z -ρ 1 ) 2 = L 2 ; (x + 1/2) 2 + (y - √ 3/2) 2 + (z -ρ 2 ) 2 = L 2 ; (x + 1/2) 2 + (y + √ 3/2) 2 + (z -ρ 3 ) 2 = L 2

Singularity Analysis

Singularities of a robotic manipulator are important feature that essentially influence its capabilities. Mathematically, a singular configuration may be defined a rank deficiency of the Jacobian describing the differential mapping from the joint space to the workspace and vice versa. Differentiating the constraint equations of the robot with respect to time leads to the velocity model: At + B q = 0 where A and B are the parallel and serial Jacobian matrices, respectively, t is the velocity of P and q defines the joint velocities. The parallel singularities occur whenever det(A) = 0 and the serial singularities occur whenever det(B) = 0 [START_REF] Chablat | Working Modes and Aspects in Fully-Parallel Manipulator[END_REF].

det(A) o = -8ρ 1 ρ 2 ρ 3 + 8ρ 1 ρ 2 z + 8ρ 1 ρ 3 y + 8ρ 2 ρ 3 x det(A) h = -8ρ 1 ρ 3 x + 8ρ 2 ρ 3 x -8ρ 1 ρ 3 + 8ρ 1 z -8ρ 2 ρ 3 + 8ρ 2 z + 16ρ 3 y det(A) t = 8ρ 1 z + 8ρ 2 z -16ρ 3 z det(A) u = 4 √ 3(3z -ρ 1 -ρ 2 -ρ 3 + ρ 3 x + ρ 2 x -2ρ 1 x) + 12ρ 3 y -12ρ 2 y (3)
Parallel and serial singularities as well as their projections in workspace and joint space are computed using a Gröbner based elimination method. This usual way for eliminating variables (see [START_REF] Cox | Using Algebraic Geometry[END_REF]) computes (the algebraic closure of) the projection of the parallel singularities in the workspace. In the same way, one can compute (the algebraic closure of) the projection of the parallel singularities in the joint space. Both are then defined as the zero set of some system of algebraic equations and we assume that the considered robots are generic enough so that both are hypersurfaces. det(A) o , det(A) h , det(A) t and det(A) u are the parallel singularities of Orthoglide, Hybridglide, Triaglide and UraneSX, respectively. Starting from the constraint equations and the determinant of the Jacobian matrix, we are able to eliminate the joint values. This elimination strategy is more efficient than a cascad-ing elimination by means of resultants which might introduce many more spurious solutions : singular points that are not projections of singular points. Due to the lack of space, other equations associated with serial singularities are not presented in this paper. Figure 2 represents the projections of parallel singularities in projection space (x, y, z). These surfaces represent the singularity associated with the eight working modes. In Table 1, a comparative study of five parameters among the family of delta like robot is presented. We have tabulated the main characteristics of the polynomials (In three variables) used for the plots (Implicit surface) : their total degree, their number of terms and the maximum bitsize of their coefficients. We have also reported the time (In seconds) for plotting the implicit surface which they define and the number of cells computed by the CAD, as well as the number of cells in the final result after gluing those that are adjacent and belongs to the same connected component. Several functions are used which involves the discriminant variety, Gröbner bases and CAD computations, computed in Maple 18 with a Intel(R) Core(TM) i7-3770 CPU @ 3.40 GHz (14 Gb RAM). As can be seen from Table 1, there exists higher values of all the parameters for the Hybridglide, among all manipulators listed, which infers that it has more complex singularities, whereas for the Triaglide all the values are least which intuits the less complicated singularities. For example, the computation times for the Hybridglyde for parallel singularities is high compared to the one for the Othoglide, even if the surface has similar characteristics. This is due to the geometry of the surface which is more difficult to decompose in the case of the Hybridglide : the Cylindrical Algebraic Decomposition is described by 1208 cylindrical cells in the case of Hybridglide while it is described by 272 cells for the Orthoglide.

x y z -2 -2 -2 2 2 2 (a) x y z -2 -2 -2 2 2 2 (b) x y z -2 -2 -2 2 2 2 (c) x y z -2 2 -2 -2 2 2 (d)

Workspace analysis

The workspace analysis allows to characterize of the workspace regions where the number of real solutions for the inverse kinematics is constant. A CAD algorithm is used to compute the workspace of the robot in the projection space (x, y, z) with some joint constraints taken in account. The three main steps involved in the analysis are [START_REF] Chablat | Workspace and joint space analysis of the 3-RPS parallel robot[END_REF][START_REF] Lazard | Solving parametric polynomial systems[END_REF][START_REF] Moroz | Sur la décomposition réelle et algébrique des systèmes dépendant de paramètres[END_REF]:

• Computation of a subset of the joint space (resp. workspace) where the number of solutions changes: the Discriminant Variety . • Description of the complementary of the discriminant variety in connected cells:

the Generic Cylindrical Algebraic Decomposition (CAD). • Connecting the cells belonging to the same connected component in the counterpart of the discriminant variety: interval comparisons.

The different shapes of workspace for the delta-like robots is shown in Fig. 3, where blue, red, yellow and green regions correspond to the one, two, four and eight number of solutions for the IKP. A comparative study is done on the workspace of the family of delta-like manipulator and the results are shown in Fig. 3. All the workspace are plotted in the rectangular box, where [START_REF] Chablat | Workspace and joint space analysis of the 3-RPS parallel robot[END_REF], so that the shapes of these workspace can be compared. From the Fig. 3 it can be intuited that the Triaglide will be good selection, if the task space is more in horizontal plane, whereas the Orthoglide is good for the three dimensional task space.

x ∈ [-2, 2], y ∈ [-2, 6] and z ∈ [-2,

Conclusions

A comparative study on the workspace of different delta-like robots gives the idea about shape of the workspace, which further plays an important role in the selection of the manipulator for the specific task or for the trajectory planning. The main characteristics associated with the singularities are tabulated in Table 1, which also gives some information about the complexity of the singularities, which is an essential factor for the singularity-free path plannings. From these data, it can be observed that the singularities associated with the Hybridglide are complicated, whereas the structure of those associated with the Triaglide is rather simple.

Fig. 1

 1 Fig. 1 Configuration plot for Orthoglide (a), Hybridglide (b), Triaglide (c) and UraneSX (d) robot

Fig. 2

 2 Fig. 2 Projection of parallel singularities of Orthoglide (a), Hybridglide (b), Triaglide (c) and UraneSX (d) in workspace. All the computations are done with joint limits.

Fig. 3

 3 Fig. 3 Workspace with joint constraints for Orthoglide (a), Hybridglide (b), Triaglide (c) and UraneSX (d) robot with the number of inverse kinematic solutions. Blue, red, yellow and green regions corresponds to the one, two, four and eight number of solutions for the IKP.

Table 1

 1 Comparison of the different parameters associated with the singularities for the robots

	Manipulators Types of	Plotting Time Degrees No.of terms Binary No. of Cells
		Singularities	(s)		Size
	Orthoglide Parallel	0037.827 18[10,10,10]	097	015 [02382,0272]
		Serial	0005.133 18[12,12,12]	062	012 [00044,0004]
	Hybridglide Parallel	5698.601 20[16,08,12]	119	017 [28012,1208]
		Serial	0007.007 18[12,12,12]	281	017 [00158,0027]
	Triaglide	Parallel	0010.625 03[00,00,03]	002	002 [00138,0004]
		Serial	0005.079 06[06,06,06]	042	007 [00077,0017]
	UraneSX	Parallel	0022.625 06[06,04,00]	015	040 [02795,0070]
		Serial	0018.391 12[12,12,12]	252	151 [00392,0142]

Acknowledgements The work presented in this paper was partially funded by the Erasmus Mundus project "India4EU II".


Iulia Dragomir
email: iulia.dragomir@irit.fr

Iulian Ober
email: iulian.ober@irit.fr

Christian Percebois
email: christian.percebois@irit.fr

Integrating Verifiable Assume/Guarantee Contracts in UML/SysML

The compositional approach based on components and driven by requirements is a common method used in the development of critical realtime embedded systems. Since the satisfaction of a requirement is subject to the composition of several components, defining abstract and partial behaviors for components with respect to the point of view of the requirement allows for a manageable design of systems. In this paper we consider such specifications in the form of contracts. A contract for a component is a pair (assumption, guarantee) where the assumption is an abstraction of the component's environment behavior and the guarantee is an abstraction of the component's behavior given that the environment behaves like the assumption. In previous work we have defined a formal contract-based theory for Timed Input/Output Automata with the aim of using it to express the semantics of UML/SysML models. In this paper we propose an extension of the UML/SysML language with a syntax and semantics for contracts and for the relations they must satisfy. Besides the important role that contracts have in design, they can also be used for the verification of requirement satisfaction and for their traceability.

Introduction

Nowadays critical real-time embedded systems grow larger in size and more complex. Their development is a challenging task and is often error-prone. A way for system designers to tackle this issue is to use a compositional approach driven by requirements. For example, process-oriented standards such as DO-178C [START_REF]Software Considerations in Airborne Systems and Equipment Certification[END_REF] highlight the need to model requirements at different levels of abstractions during development and to ensure their traceability at each design iteration step.

However, requirements are often difficult to be mapped to components: several components combine together to satisfy a requirement and a component may be involved in the satisfaction of several requirements. In order to achieve provably correct compositional design, one needs a way to abstractly specify how a particular component K participates in fulfilling a requirement ϕ. Such a specification can take the form of a contract: a pair (assumption, guarantee) where the assumption is an abstraction of K's environment behavior and the guarantee is an abstraction of K's behavior given that the environment behaves according to the assumption. Such a contract can then be used to model the point of view of the component with respect to the requirement ϕ. Contracts for reactive and real-time components have received a lot of attention from the research community recently, as discussed in §5.

Besides the important role contracts can play in system design, they can also be used as basic blocks for compositional verification of requirement satisfaction. In [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF] we have introduced a contract-based theory for compositional verification of systems of communicating Timed Input/Output Automata (TIOA) with the intention to use it as underlying semantics for contract-based UML/SysML [START_REF]Object Management Group[END_REF][START_REF]Object Management Group[END_REF] modeling and verification.

This paper complements the theory from [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF] by extending UML/SysML with the language elements needed for modeling contracts and their relations. The paper is structured as follows: in §2 we summarize the contract-based reasoning theory we have defined and we present the OMEGA UML/SysML profile [START_REF] Conquet | Formal Model Driven Engineering for Space Onboard Software[END_REF] on which we want to apply the formal theory. In §3 we propose a meta-model for the contract theory and a set of constraints and well-formedness rules needed to make the system model verifiable with contracts. Then, an instantiation of the meta-model for the OMEGA UML/SysML profile is discussed. §4 presents the application of our approach to an industry-grade system model, the ATV SGS case study previously described without contracts in [START_REF] Dragomir | A case study in formal system engineering with SysML[END_REF], before concluding.

Background

2.1 Timed Input/Output Automata Many mathematical formalisms have been proposed in the literature for modeling communicating timed reactive components. Our work is based on a variant of Timed Input/Output Automata of [START_REF] Kaynar | The Theory of Timed I/O Automata -Second Edition[END_REF] since it is one of the most general formalisms, thoroughly defined and for which several interesting compositionality results are already available. A TIOA specifies a state space and a set of admitted timed behaviors for a component. The parallel composition of TIOAs (denoted in the following) is based on synchronization of corresponding inputs/outputs and the interleaving of other actions. The main differences between our variant and that of [START_REF] Kaynar | The Theory of Timed I/O Automata -Second Edition[END_REF] are that (1) matched input-outputs resulting from a composition of two automata result in a visible I/O action which is not involved in synchronizations thereafter (in [START_REF] Kaynar | The Theory of Timed I/O Automata -Second Edition[END_REF] the result is an output, which is not consistent with the semantics of signals in UML or SysML) and (2) continuous variables are restricted to Alur-Dill-style linear clocks [START_REF] Alur | A theory of timed automata[END_REF] in order to allow the use of known symbolic simulation and verification methods (unavailable for the more general model of [START_REF] Kaynar | The Theory of Timed I/O Automata -Second Edition[END_REF]). For the full definition of the formal model, the reader is referred to [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF].

A UML/SysML Profile Based on Timed Input/Output Automata

In previous work [START_REF] Bozga | The IF Toolset[END_REF] we have considered the high-level modeling of embedded realtime systems in UML/SysML with a semantics provided in terms of TIOA. The result of this work is a semantic profile for UML and SysML called OMEGA and a set of tools for simulation and model-checking, the IFx-OMEGA toolset 1 .

In OMEGA, the architecture of a system is expressed in the usual way in UML (class diagrams) and in SysML (block definition diagrams and internal block diagrams). Classes/blocks may use most of their features: properties (attributes and parts), signals receptions, interconnection elements (port, connector, interface) and relations (association, composition and generalization). The hierarchical architecture of components (and systems) is specified through composite structures.

The behavior of atomic components is modeled by state machines with standard UML actions on transitions. The operational semantics of each component instance is a timed input/output automaton. The TIOAs corresponding to components are composed in parallel and communicate by asynchronous signal exchanges. This imposes that all communications between objects/block instances are defined as signal outputs and receptions that are transferred via ports and connectors. Ports need to be typed with interfaces that contain the list of signals transferred to or from the component's environment. With respect to the actions described on transitions, the semantics supports a significant part of the UML action meta-model, including signal output, assignment, expression valuation and control flow structuring statements. The translation between OMEGA modeling concepts and the underlying TIOA semantics is automatic in IFx-OMEGA.

Figure 1: Contract-based reasoning for a subsystem with three components [START_REF] Quinton | Design, vérification et implémentation de systèmes à composants[END_REF].

The temporal elements of TIOA, such as clocks, clock actions and timed guards actions have corresponding language constructs in the OMEGA profile. The elapsing of time is constrained by transition urgency stereotypes inspired from timed automata with urgency [START_REF] Bornot | An algebraic framework for urgency[END_REF]: time delay is blocked if one of the active transitions in the current state is stereotyped «eager », it is upper-bounded if an active transition is stereotyped «delayable» and is unbounded otherwise (i.e., if all active transitions are stereotyped «lazy», which is by default).

The profile also proposes mechanisms for formalizing requirements, in particular in the form of real-time safety properties described by observer classes (identified by a stereotype «observer »). The state machine of these classes uses special primitives for monitoring the system state and events and gives verdicts about the (non-)satisfaction of a property by using labels (e.g., stereotype «error ») on states.

For a more complete description of the UML/SysML component model used in OMEGA and of the mapping of notions, the reader is referred to [START_REF] Ober | OMEGA2: A New Version of the Profile and the Tools[END_REF].

A Theory of Contracts for Timed Input/Output Automata

In [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF] we defined a theory for modeling and reasoning with contracts for TIOA. The theory is an extension of a meta-theory defined in [START_REF] Quinton | Contract-based verification of hierarchical systems of components[END_REF][START_REF] Quinton | Design, vérification et implémentation de systèmes à composants[END_REF]. A contract for a TIOA component K is a pair of TIOAs that model the assumption (A) and the guarantee (G). The satisfaction of the contract is defined formally by a relation based on trace inclusion between K A and G A modulo the set of actions that are of interest for the contract. The theory also provides the necessary mechanisms for compositional reasoning with contracts, explained in the following.

Consider that the objective is to prove that a system S composed of several components K 1 , K 2 , • • • , K n satisfies a property ϕ (see Fig. 1) under a certain hypothesis A on the behavior of its environment. The method consists in defining a more abstract specification G of the system such that A G satisfies ϕ (conformance step in Fig. 1). Once this step performed, one will be inclined to verify that the system S satisfies the contract C = (A, G). However, it is often impossible to verify directly (i.e. by checking trace set inclusion) that the composite system S satisfies the contract because of the combinatorial explosion of the state space. To avoid this problem, the method defined in [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF][START_REF] Quinton | Design, vérification et implémentation de systèmes à composants[END_REF] uses a decomposition of the proof in independent steps based on the definition of a set of individual contracts C 1 , C 2 ... for Figure 2: UML meta-model extended with contracts. the components K 1 , K 2 ..., which, when put together, ensure the global contract C. We say that {C 1 , C 2 , • • • } dominates C (dominance step in Fig. 1). The theory in [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF][START_REF] Quinton | Design, vérification et implémentation de systèmes à composants[END_REF] provides a set of sufficient conditions for dominance which can be checked independently with lesser combinatorial complexity. In addition to the conditions for dominance, one also has to check that each component K i satisfies the contract C i (satisfaction step in Fig. 1).

Extending UML and SysML for Modeling Contracts

In this section we present the UML and SysML extensions that we propose in order to support modeling and reasoning with contracts. In §3.1 we describe a domain meta-model of the contract-related concepts which extends the UML meta-model. The adaptation for SysML is discussed briefly at the end, as is relatively straightforward. We then discuss in §3.2 the constraints and well-formedness rules imposed on the key notions in order to make models with contracts compliant with the theory from [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF] and thus verifiable. Finally, in §3.3 we discuss the mapping of the metamodel concepts as a UML/SysML profile, using the standard extension mechanisms (stereotypes).

A Meta-Model for Contracts in UML

Within the contract theory we have presented there are two categories of concepts:

(1) those related to contract modeling represented in the upper part of the metamodel given in Fig. 2 and (2) those related to modeling relations between contracts, used for example in verification, represented in the lower part of Fig. 2.

The requirement ϕ that the component model has to satisfy is represented by the meta-class SafetyProperty. This meta-class is left unspecified at this point since different formalisms could be used to model a requirement, such as temporal logics, automata-based languages, etc. In §3.3 we instantiate the meta-model in the OMEGA profile, which uses an observer for modeling a SafetyProperty.

The assumption/guarantee of a contract is modeled by the corresponding metaclass Assumption/Guarantee type of Class. Each of these two elements is modeled by a class that has a behavior expressed by a state machine and communicates (only) through ports. Assumptions/guarantees are thus described in the same language as the system components. In order to have a clear semantics in terms of TIOA, they are subject to few restrictions: they should not be involved in associations, generalizations, realizations (except the interface realizations demanded by ports) and dependencies. An assumption or guarantee may define a composite sub-structure. Such an example is provided in the case study of §4.

A contract is represented by the meta-class Contract as a composite structure, containing exactly one assumption and one guarantee (i.e. any other properties are forbidden) and does not exhibit any behavior. The only relations a contract may be involved in are those that represent the verification relations used in our theory, as described below.

If a contract serves in the conformance step in the methodology depicted in Fig. 1, this is modeled using a Conformance relation (a kind of Dependency) between the Contract and the corresponding SafetyProperty. One can use the same contract for several safety properties.

The dominance relation is represented by the meta-class Refinement type of Dependency. One contract is refined by a set of contracts. Note that this is possible since UML defines Dependency from n clients to n suppliers. To ensure that no cycles may be modeled, the following constraint is imposed: the target of a Refinement is not a member of the source set.

Finally, the relation between a component and the contract that it must satisfy is represented by two relations: one at the level of the type of the component and one at the level of the instance (the part which participates in a composite structure where the contract is relevant and which is modeled by the meta-class Property). On the level of the type, an Implementation relation (a kind of Dependency) between a class and a contract models the fact that the class satisfies the contract. One class can satisfy several distinct contracts. On the level of instances, a ContractUse relation (also a kind of Dependency) between a Property, which is part of a composite structure, and a Contract models the fact that the contract is used for verification within the context of that composite structure. A Property may use a contract if and only if its class implements that contract.

All constraints have been formalized with OCL [START_REF]Object Constraint Language (OCL)[END_REF]. Their code can be found in an extended version of this paper [START_REF] Dragomir | Integrating verifiable Assume/Guarantee contracts in UML/SysML[END_REF]. This meta-model can be easily adapted for the SysML language by defining the meta-classes Contract, Assumption and Guarantee as extensions of the stereotype block applied on the meta-class Class from UML4SysML package. Similarly, the meta-class Dependency belongs to the same package.

Well-Formedness Rules for Verifiable Contracts

In order to be able to apply the contract-based verification theory from [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF] we need to make sure that the hypotheses and constraints imposed by the formal framework are satisfied by the system model. In the following we formalize these constraints at the meta-model level by a set of well-formedness rules.

Within the formal framework, a contract is modeled by a pair (A, G) of TIOA such that the set of inputs/outputs of G is a subset of the set of inputs/outputs of the component implementing the contract and the composition of A and G is a closed system. To ensure this, the set of ports of a Guarantee must correspond to a subset of the set of ports of the component for which the guarantee is defined. The correspondence is based on the port name and the port type and direction must coincide. We consider that when a port is present in the guarantee, all the corresponding signal receptions defined by the port type are handled in the guarantee. Moreover, as the composition between an Assumption and a Guarantee must be closed, every port of the Guarantee must have a corresponding conjugated port on the side of the Assumption, with the same type and reversed direction.

The dominance relation is also subject to the refinement of provided/required requests. This rule is also expressed with respect to ports: a port of the guarantee which is the target of the refinement must be matched (by name and type) by a port of one of the refining guarantees and must not be matched by a corresponding conjugated port (i.e. with reversed directionality) of another from the refining guarantees.

The theory from [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF] also induces some constraints on the state machines of assumptions and guarantees. In particular, the behavior of an assumption or guarantee should not impose constraints on time progress. This is realized on the UML/SysML level ensuring that all transitions in these state machines are stereotyped «lazy» and that there is at most one output action on any transition.

Furthermore, for a model with contracts to be used in compositional verification according to the methodology described in §2, the model must describe a unique and complete proof tree: all implemented contracts are used within a context and for all SafetyProperty there is a contract conforming to it. The rules described above have been formalized in OCL and can be found in [START_REF] Dragomir | Integrating verifiable Assume/Guarantee contracts in UML/SysML[END_REF].

Instantiating the Meta-Model in the OMEGA Profile

In order to use contracts in a standard UML or SysML model, one needs to capture the information from the meta-model described in the previous section in the form of standard extensions, namely using stereotypes. Since all the new concepts introduced in the meta-model inherit from an existing meta-class (either Class or Dependency), we choose to represent them as stereotypes of these base meta-classes. Thus, we use for the meta-class Class the stereotypes «contract», «assumption», «guarantee», and «observer » (which already exists in OMEGA and is reused for representing SafetyProperty). For contract relations, the stereotypes of Dependency that correspond to the meta-model elements are «contractImplementation», «con-tractUse», «contractRefinement» and «contractConformance».

As explained before, from the semantic point of view, Contracts are not handled in the same way as usual classes/blocks as they are not considered executable elements of the system. Contracts are only used by the verification tools to check the validity of the conformance, dominance and satisfaction relations. The choice to describe contracts as classes is due to the fact that the syntax of contracts reuses much of the standard syntax of classes: we need to represent a contract as a closed composite structure with ports and links for the communication between the assumption and the guarantee and between the assumption and the component which uses the contract, interfaces to statically type the sets of inputs and outputs of each sub-component and behavior in the form of state machines.

The ATV Solar Generation System Case Study

The concepts and the reasoning method previously described have been applied on a case study, an industrial-grade system model of a subsystem of the Automated Transfer Vehicle (ATV). The ATV, developed by Astrium Space Transportation for the European Space Agency, is a spacecraft put into orbit by the European heavy launcher Ariane-5 with the aim of supplying the International Space Station. This case study consists of the Solar Wing Generation System (SGS) [START_REF] Dragomir | A case study in formal system engineering with SysML[END_REF] responsible for the deployment and management of the solar wings of the vehicle. The SysML model used in the following, provided by Astrium Space Transportation, was obtained by reverse engineering the actual SGS system for the purpose of this study.

The SGS system model 2 illustrated in Fig. 3 summarizes the three main components involved in the case study: the mission and vehicle management (MVM ) part that initiates SGS wing deployment, the SOFTWARE part of the SGS that based on requests received from the MVM executes the corresponding automated procedures and the HARDWARE part that models the four physical wings. The communication between components is realized via asynchronous signals transported through ports and connectors. Due to the large number of ports (661) and connectors (504), Fig. 3 presents a simplified architectural view of SGS and only shows a link between two parts where several connectors and ports are involved in the actual model. Let us mention that the verification steps presented below have been performed on the initial system model.

Under the hypothesis that at most one hardware failure may occur during a run, which is embedded in the HARDWARE model, the main goal of the case study is to verify the following property ϕ: after 10 minutes from system start-up, all four wings are deployed.

Due to the size and the complexity of the model, applying model-checking directly leads to combinatorial explosion and the verification of ϕ does not finish. We explain in the following how the property ϕ was verified using the contract-based reasoning methodology. We start by modeling the property ϕ. This implies identifying what the observer corresponding to the safety property ϕ must monitor. In our case, the block phi must observe the answer that each wing provides with respect to its status (deployed or not deployed) when interrogated by the software. So, the property ϕ expressed with respect to wing behavior must be satisfied by the HARD-WARE block instance that contains them. With regard to Fig. 1, the HARDWARE is the subsystem S and WINGi, i = 1, 4, are the components K i . The environment of the subsystem is given by the parts with which it communicates: bidirectional communication is directly established between SOFTWARE and HARDWARE, while SOFTWARE depends on the behavior of MVM. Thus, the environment E of Fig. 1 is represented here by the composition of MVM and SOFTWARE.

Next, we provide a contract C = (A, G) such that it conforms to ϕ. In order to comply to the contract methodology, C is implemented by HARDWARE 's type and it is used by this part within the proof tree. We use as assumption A the concrete environment of HARDWARE, the composition between MVM and SOFTWARE itself, which thus satisfies by construction the mirror contract C -1 = (G, A). Keeping this composition as assumption is not problematic since its state space has a manageable size. As guarantee G we use the following abstraction derived (manually) from the individual behavior of wings: for each wing status interrogation, the target wing answers either as not deployed for at most 400 seconds or as deployed after at least 130 seconds. In order to ensure that the contract defines a closed system, since MVM SOFTWARE sends all possible requests to HARDWARE, we equip G with all ports defined by HARDWARE and we enrich the behavior of G to ignore all other requests. Then G and HARDWARE have the same set of ports: no refinement of requests is performed and the corresponding rule is satisfied.

The third step consists in modeling a set of contracts {C 1 , C 2 , C 3 , C 4 } that refine C and proving that each contract

C i = (A i , G i ) is implemented by WINGi 's type, i = 1, 4.
The environment for WINGi is given by the environment of the subsystem HARDWARE and all WINGj, j = i. We use the following abstraction WAj for WINGj : the wing is either not deployed for at most 400 seconds or deployed from at least 130 seconds while all other received requests are consumed. The assumption A i is the parallel composition of MVM, SOFTWARE and WAj, j = i. The guarantee G i 2 The case study is represented using the notation conventions of IBM Rhapsody: http: //www-03.ibm.com/software/products/us/en/ratirhap/. is the projection of G on WINGi. Again the refinement of requests is not considered and since the ports of HARDWARE with respect to WINGi are identical to the ports of WINGi, both rules on port set inclusion for implementation and refinement are satisfied. Moreover, each C i defines a closed system.

After this step the proof obligation tree is complete. The verification involves 10 intermediate steps: 4 for verifying that each wing satisfies its contract, 5 sufficient conditions for dominance between {C 1 , C 2 , C 3 , C 4 } and C and one for proving that A G ϕ. Each verification step is performed by the OMEGA-IFx model checker in a few hours. The overall effort of the building the contracts and performing the verification steps was of about 5 person*days. Even though the translation and verification algorithm are automated, two intermediate modeling steps -connecting assumptions to components via links and transforming the state machine of guarantees into observer (timed trace inclusion is verified using observers) -remain manual. The automation of these steps is currently under development. For further details on contract-based verification of the SGS case study, the reader is referred to [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF].

Related Work

Modeling and verifying contracts for components is a long line of research, whose origins date back to Hoare logic [START_REF] Hoare | An Axiomatic Basis for Computer Programming[END_REF]. Syntactical and behavioral contracts, as classified in [START_REF] Beugnard | Making Components Contract Aware[END_REF], have been explored for specifying composition constraints and pre/post conditions for operations and also for modeling transformation of models and execution semantics. Contracts as a language construct have emerged with the Eiffel programming language [START_REF] Meyer | Applying Design by Contract[END_REF] and have since been explored for various programming and specification models. In this section we concentrate on work aiming to introduce contracts in high-level modeling languages. For a discussion of more theoretical works on contracts and contract-based verification the reader is referred to [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF].

Weis et al. [START_REF] Weis | A UML Meta-model for Contract Aware Components[END_REF] propose to model a contract for a component in UML by an interface and to specify its role: it can be either a required contract on which the component depends or a provided contract that is realized by the component. Syntactically, this representation of contracts is similar to ours: we also make the distinction between the required behavior of the environment and the provided behavior of the component by taking into account the assumption over the environment. However, our contracts are richer since they model a behavior that can be used for component validation, while the contracts of [START_REF] Weis | A UML Meta-model for Contract Aware Components[END_REF] can be used only for composability checking during the development phases.

The Kmelia component model [START_REF] Messabihi | Multilevel Contracts for Trusted Components[END_REF][START_REF] André | Vérification de contrats logiciels a l'aide de transformations de modeles[END_REF], based on the work described above, provides means to verify the functional correctness of behavioral contracts for services: the behavior of an operation is modeled as a Labeled Transition System and formal verification can be realized within different tools via model transformation. Their meta-model defines for a contract the source implementing it as an aggregated element (operation or interface) and models explicitly the contract satisfaction results. But this formalism does not describe how the order in which services are called by and from a component can be verified, order that can be seen similar to our state machines from components. Furthermore, it does not provide a connection to high-level modeling languages as UML/SysML.

Contracts modeled as pre/post conditions are used in [START_REF] Cariou | OCL contracts for the verification of model transformations[END_REF] for the verification of model transformation: the assumption is represented by an OCL constraint on the source model and the guarantee is an OCL constraint on the target model. Moreover, with respect to the syntax of contracts the two approaches are different: while the one we describe considers contracts only for components, [START_REF] Cariou | OCL contracts for the verification of model transformations[END_REF] models contracts for all model elements. The same contracts are used in [START_REF] Cariou | Contracts for model execution verification[END_REF] to model the execution semantics of UML elements which is seen as a particular case of model transformation.

To the best of our knowledge, this study is the first to consider behavioral contracts at the component level in UML/SysML and to provide verification relations for property satisfaction by contract-based reasoning. The meta-model we propose is generic enough to represent all the other meta-models previously described, excepting the verification results extension that is based on the dynamical execution of the model.

Conclusion

Based on a theory of contracts and on a methodological approach for reasoning with contracts introduced in previous work [START_REF] Dragomir | Safety Contracts for Timed Reactive Components in SysML[END_REF], we have proposed an extension of UML/SysML allowing to model contracts and use them for compositional verification of requirements. The extension is defined as a meta-model, enriched with constraints and well-formedness rules to make contracts verifiable. We have instantiated the extension within the OMEGA UML/SysML profile to make it usable with standard model editors. The verification method is supported by the OMEGA-IFx toolset and the approach was validated on an industrial-grade system model.

Although an automatic model transformation from OMEGA system models to the input language of the IFx Toolset is already available, some of the steps for generating the intermediate contract-based verification models remain manual. Future work consists in automating all the intermediate model generation steps and in adding functionality for managing the proof obligations and results and for enforcing the rigorous verification methodology described in §2.3.

The method described before does not explicitly prescribe how to derive contracts for the whole system and for its components. In the case study described in §4, this task was relatively straightforward: since we make no additional assumption (A) about the environment, G is roughly the same as ϕ, and the component guarantees are a projection of the desired global guarantee. There may be cases where the definition of contracts is less obvious and the overhead is significant, and previous attempts to introduce contracts in programming have not enjoyed a widespread success due to this kind of overhead. Nevertheless, we believe that the case for contracts in early phases of system engineering is different than the case for software programming, and the overhead should be acceptable for certain critical systems. Further work is needed in order to lower the overhead by finding methods or methodological guidelines for deriving intermediate contracts from the properties one is trying to prove.

Figure 3 :

 3 Figure 3: Architecture of the SGS system including contracts (simplified view).

http://www.irit.fr/ifx


Manuel Marin

David Defour

Linear circuit analysis based on

come

Linear circuit analysis based on parallel asynchronous fixed-point method Manuel Marin, Student Member, IEEE, David Defour, and Federico Milano, Senior Member, IEEE Abstract-Time-domain circuit analysis and simulation is commonly performed through intrinsically serial methods, such as those implemented in the popular software package SPICE. However, the practical utilization of these methods is limited to circuits of a certain size, as for larger problems the simulation time becomes prohibitive. Parallelization of such serial routines has been proposed as the main alternative to accelerate the analysis and overcome the problem, even if the speed-up that can be achieved by this strategy is bounded, according to the Amdahl's law. Currently, there is a lack of intrinsically parallel methods which would allow to detach the efficacy of the solution from the problem size. In this article, we develop a new theoretical approach to circuit analysis from an intrinsically parallel point of view. We propose a fixed-point method and determine its convergence condition according to the theory of asynchronous iterations. We also perform a series of tests that show that our method is faster in most cases than those based on the traditional intrinsically serial approach. In particular, we obtain an empirical proof that our approach is independent of the problem size, offering great opportunities for scalability.

Index Terms-Linear circuits, time-domain analysis, fixedpoint arithmetic, parallel algorithms, iterative method.

I. INTRODUCTION

C IRCUIT simulation has become an extremely important part in integrated circuit evaluation and design and a continuously growing field [START_REF] Brambilla | The simulation errors introduced by the spice transient analysis[END_REF], [START_REF] Brambilla | Envelope following method for the transient analysis of electrical circuits[END_REF], especially after the introduction of the SPICE simulator [START_REF] Nagel | SPICE: Simulation program with integrated circuit emphasis[END_REF]. Meanwhile, as the number of components in modern applications continues to expand, there is a constant necessity to analyze larger and larger circuits, eventually without increasing the simulation time [START_REF] Elwan | Low-voltage low-power cmos current conveyors[END_REF], [START_REF] Awwad | On modeling of parallel repeater-insertion methodologies for soc interconnects[END_REF]. Most attempts to achieve such goal have consisted of taking one of the existing methods, which are intrinsically serial, and parallelize part of it in order to obtain a speed-up [START_REF] Gulati | Fast circuit simulation on graphics processing units[END_REF]- [START_REF] Peng | Parallel transistor level circuit simulation using domain decomposition methods[END_REF]. Nevertheless, in doing this, there is always a fraction of the algorithm that remains sequential, which according to the Amdahl's Law determines that the speedup to be achieved is fairly limited [START_REF] Amdahl | Validity of the single processor approach to achieving large scale computing capabilities[END_REF]. This would not be the case of an algorithm that is intrinsically parallel; however, the task of developing such algorithm requires to reformulate the way in which we analyze circuits. In this article, we present a novel approach to accelerate circuit analysis and simulation which instead of parallelizing existing, intrinsically serial routines, proposes a new, intrinsically parallel method allowing to M. Marin is with Université de Perpignan Via Domitia, DALI and Université Montpellier 2, LIRMM, France, and also with the School of Electrical, Electronic and Communications Engineering of the University College Dublin, Dublin, Ireland (e-mail: manuel.marin@univ-perp.fr). D. Defour is with Université de Perpignan Via Domitia, DALI and Université Montpellier 2, LIRMM, France (e-mail: david.defour@univ-perp.fr).

F. Milano is with the School of Electrical, Electronic and Communications Engineering of the University College Dublin, Dublin, Ireland (e-mail: federico.milano@ucd.ie).

detach the complexity of the solution from the problem size. The proposed model attempts to change circuit simulation paradigm and is inspired from the circuit component dynamics that can be observed in real life.

The transient response of a circuit is typically obtained as a series of system states computed for several time instants. Computing one state involves the resolution of the circuit equations evaluated at that point. In SPICE, this is done in two sequential steps, load and solve. During the load step, the circuit equations are obtained from evaluating component models; during the solve phase, a linear system is solved. As component models are generally non-linear, the resolution of the system equations usually involves the repetition of the two steps in an iterative procedure, such as the Newton-Raphson method. In some applications, the load phase takes up to 75% of the total time, although this ratio tends to decrease with the circuit size.

The load phase is easy parallelizable: it is done by distributing the evaluation of several independent components among independent processors. Some successful examples of this, involving Field-Programmable Gate Arrays (FPGA) and Graphic Processing Units (GPU) are presented in [START_REF] Gulati | Fast circuit simulation on graphics processing units[END_REF], [START_REF] Kapre | Performance comparison of single-precision spice model-evaluation on fpga, gpu, cell, and multi-core processors[END_REF]. The solve phase, on the other hand, requires a more delicate approach. The original version of SPICE uses sparse matrix routines to solve the linear system, and PARASPICE is one of the first attempts of parallelizing such routines [START_REF] Yang | Paraspice: a parallel circuit simulator for shared-memory multiprocessors[END_REF]. More recently, in [START_REF] Kapre | Parallelizing sparse matrix solve for spice circuit simulation using fpgas[END_REF], the same approach is revisited using the KLU algorithm [START_REF] Davis | Algorithm 907: Klu, a direct sparse solver for circuit simulation problems[END_REF] and FPGA. Other researchers have proposed relaxation techniques (e.g., Gauss-Seidel, Gauss-Jacobi) as an alternative, more parallel-friendly approach to matrix factorization. A survey of those techniques, including timing simulation, iterative timing simulation and waveform relaxation is done in [START_REF] Newton | Relaxation-based electrical simulation[END_REF]. In [START_REF] Peng | Parallel transistor level circuit simulation using domain decomposition methods[END_REF], a novel approach consisting of decomposing the domain and distributing it into several parallel processors is presented. However, at the end, sequentiality persists in all these models, thus, limiting the maximum speed-up according to Amdahl's Law.

The approach that we develop in this article attempts to avoid sequentiality by mimicking the actual behaviour of an electrical circuit, where each component evolves in a continuous cycle, permanently reevaluating its state regardless of external factors. To efficiently model such a process, we apply the concept of team algorithms, developed in [START_REF] Talukdar | Designing Algorithms and Assignments for Distributed Processing[END_REF]. In this concept, the value of a critical unknown is computed through running several different algorithms and forming convex combinations of the results; the obtained value is used as starting guess in a new cycle, until all the algorithms reach an agreement on the unknown's value. An analog idea is presented in [START_REF] Tsitsiklis | Problems in decentralized decision making and computation[END_REF] under the name of agreement algorithm.

Similarly, in the real operation of an electrical circuit, we have a number of electrical components connected through a common node; the components impose their dynamics on the node, affecting and modifying the node's voltage, until it becomes stable after a certain number of interactions. We can argue that this kind of physical phenomena that occur in real life are intrinsically parallel.

This paper intends to bring the following main contributions to the existing literature:

1) It proposes a novel method for time-domain transient circuit simulation, which is based on an intrinsically parallel iteration and so is unaffected by the Amdahl's argument. A proof of convergence of the iterative method is given for the case of linear component models.

2) It presents an empirical comparison of the proposed method with traditional approaches based on matrix factorization. In particular, it is shown that the number of iterations performed by our method is independent from the problem size, allowing performance to scale with the amount of computing units available. The remaining of the article is organized as follows. In Section II, we offer a discussion on several alternatives to implement an iterative process on a parallel computer, which serves as theoretical background for the following sections. In Section III, we present our model of intrinsically parallel circuit simulator and our convergence result. Implementation details and the results of a series of tests, along with comparisons with existing circuit simulators based on the traditional, serial approach, are given in Section IV. Finally, Section V concludes the paper.

II. BACKGROUND

Usually, iterative processes are implemented on a parallel computer by assigning each component of the system state vector to an independent processor. From that point on, different synchronization schemes lead to different types of iteration. The first one that we want to revisit is the synchronous iteration, defined as follows.

Definition 1 (Synchronous iteration). Let x ∈ R n and f : R n → R n . Given an initial value x 0 , the series defined by:

x k+1 = f (x k ), k ∈ N, (1)
is termed a synchronous iteration associated to F (•).

In this model the vector is updated as a block, i.e., all the components are updated before the process moves to a next iteration. In order to implement it, a synchronization step is needed at the end of each cycle. In this synchronization step, processors that have finished their update wait for the others to catch up. The idle time will depend on a series of factors, e.g., load balancing, processor frequency, communication network and memory access patterns [START_REF] Frommer | On asynchronous iterations[END_REF].

An alternative to the synchronous model is the so-called asynchronous iteration, in which individual processors are always allowed to work, independently of the state of others. Typically, one processor will update the component assigned to it using the most recent data available, then communicate its action to other processors and restart the cycle, without waiting. In consequence, only certain vector components are updated at each asynchronous iteration, and some of the components used in these updates have values that do not correspond to the very last iteration. In order to specify the model, we need to introduce the update function and the delay function.

The update function, noted U(•), receives the iteration counter k ∈ N, and returns a subset of {1, . . . , n} indicating the list of processors that will update their components during iteration k. The delay function, noted d(•), receives the indices i, j ∈ {1, . . . , n} and the iteration counter k ∈ N, and returns the delay of processor j with respect to processor i at iteration k.

Definition 2 (Asynchronous iteration). Let x ∈ R n and f : R n → R n . For k ∈ N, i, j ∈ {1, . . . , n}, let U(k) ⊆ {1, . . . , n} and d(i, j, k) ∈ N 0 , such that:

d(i, j, k) ≥ 0, ∀i, j, k, (2a) lim k→∞ d(i, j, k) < ∞, ∀i, j, (2b)
|{k : i ∈ U(k)}| = ∞, ∀i. (2c)
Given an initial value x 0 , the series defined by:

x k+1 i = f i (x k-d(i,1,k) 1 , . . . , x k-d(i,n,k) n ) if i ∈ U(k), x k i if i ∈ U(k), (3)
is termed an asynchronous iteration associated to f (•), with update function U(•) and delay function d(•).

The assumption in (2a) states that only values computed in previous iterations are used in any update. The one in (2b) states that newer values of the components are eventually used. Finally, the assumption in (2c) states that no component ceases to be updated during the course of the iteration. In the case of a shared memory machine, the two latter assumptions become equivalent, as there are no delay associated to the communications. Indeed, as soon as a component is updated, its new value becomes visible and available to all the processors. Also; note that a synchronous iteration is a special case of asynchronous iteration with U(k) = {1, . . . , n} and d(i, j, k) = 0, for all i, j, k.

Regarding convergence of asynchronous iterations, several results have been published. In particular, when the application function f (•) is linear we have the following theorem by Chazan and Miranker [START_REF] Chazan | Chaotic relaxation[END_REF].

Theorem 1 (Sufficient condition for convergence in the asynchronous case [START_REF] Chazan | Chaotic relaxation[END_REF]). Let x ∈ R n , and f : R n → R n be a linear application, i.e.,

f (x) = Lx + b, L ∈ R n×n , b ∈ R n . ( 4
)
Let |L| denote the matrix of absolute values of the entries of L, and ρ(|L|) its spectral radius (see Appendix A for a definition of spectral radius). If

ρ(|L|) < 1, (5)
then the asynchronous iteration

x k , k ∈ N associated to f converges to x * , the unique fixed point of f , regardless of the selection of U(•), d(•) and x 0 . Furthermore, if ρ(|L|) ≥ 1,
then there exists U(•), d(•) and x 0 such that x k does not converge to x * .

In some cases, if the assumption in ( 5) is not met, one can still build a series that converges to a fixed point x * . However this requires additional assumptions on the update and shift functions, leading to a new kind of iteration known as partially asynchronous.

Definition 3 (Partially asynchronous iteration). Consider Definition 2 of an asynchronous iteration. Replace the assumptions in (2b) and (2c) by the following:

∃ d ∈ N : d(i, j, k) ≤ d, ∀i, j, k, (6a)
∃s ∈ N : i ∈ s s=1 U(k + s), ∀i, k. (6b) d(i, i, k) = 0, ∀i, k, (6c)
The series x k is now termed a partially asynchronous iteration.

The assumption in (6a) establishes that not only newer values of the components are eventually used, but each of these values is used before d iterations have passed from their calculation. The assumption in (6b), in turn, states that each component is updated at least once in every s consecutive iterations. These assumptions require to introduce a synchronization step every once in a while (for example, every d or s iterations). In the case of a shared memory machine, only one of the two is needed and d = s. The assumption in (6c) establishes that every component is updated using its last calculated value. In practical terms this is equivalent to having each component assigned to only one processor. Now we can enunciate the following result, regarding convergence of partially asynchronous iterations.

Theorem 2 (Sufficient condition for convergence in the partially asynchronous case [START_REF] Tseng | Partially asynchronous, parallel algorithms for network flow and other problems[END_REF]). Let x ∈ R n , and f : R n → R n be a linear application, i.e.,

f (x) = Lx + b, L ∈ R n×n , b ∈ R n . (7)
Also, let g : R n → R n , defined by:

g(x) = (1 -α)x + α(Lx + b), (8)
where α ∈ R and 0 < α < 1.

If L = (l ij ) is irreducible (see Appendix A for a definition of irreducibility) and n j=1 |l ij | ≤ 1, ∀i, (9)
then the partially asynchronous iteration x k , k ∈ N, associated to g(•) converges to x * , fixed point of f (•).

In [START_REF] Lubachevsky | A chaotic asynchronous algorithm for computing the fixed point of a nonnegative matrix of unit spectral radius[END_REF], Lubachevsky and Mitra presented a special case in which convergence occurs for α = 1, at a linear (geometric) rate. They also showed that the average rate of convergence per iteration in the long-term, is low-bounded by an expression which combines the problem size n, the measurements of asynchronism d and s, the entries of matrix L and the fixed point x * , as follows:

(1 -σ r ) 1/r , (10)
where

r = 1 + d + (n -1)(s + d), (11a)
σ = min i,j + x * i l ij x * j , (11b)
and min + (•) refers to the minimum of the positive elements. This means that the error cannot be reduced by a factor greater than (1σ r ) 1/r in any average iteration. Note that as σ is always lower than 1, for increasing values of n, d and s, the rate of convergence approaches 1, i.e., the convergence becomes sub-linear. Assume that we are in the synchronous case, in which d = s = 0. Then from equation (11a) we have r = 1, so the expression in (10) becomes simply

1 -σ. ( 12
)
In other words the rate of convergence is unaffected by the size of the problem, n. Now assume that we are in a shared memory environment, that only allows for the slightest desynchronization (1 iteration) between processors, i.e., d = s = 1. In this case we have r = 2n. Let's see how this affects the convergence rate: Figure 1 shows the bound on the convergence rate (lower is better) as a function of the problem size, for different values of the parameter σ. We observe that for very small problems (n < 2) the convergence is already sub-linear. Only for 'very good' values of σ (e.g., 0.9) the problem can grow, although very slightly, without the convergence rate being affected so much. This may be a strong reason to prefer the synchronous approach over the partially asynchronous. However, depending on the implementation issues, the partially asynchronous strategy may become advantageous in certain cases, e.g., when the costs of synchronization are relatively high, or when the number of iterations needed to converge is relatively low.

III. MODEL DESCRIPTION

Following the precedent discussion on parallel iterations we now introduce our model of an intrinsically parallel circuit simulator, through a top-down approach. Consider Table I as a synthesis of the variables and parameters of the model, which we introduce gradually. The general layout of a discrete-time circuit simulator is given in Algorithm 1. Instruction 5 is the crucial step, in which the system state is updated. Here is where our approach differs from SPICE and such classic techniques, in that the method applied to compute the state consists of an intrinsically parallel iteration. We will now describe this method and the structures that support it.

Algorithm 1 Discrete-time simulator.

Input: Circuit specification including initial state v 0 , time step length ∆t and total time T .

Output: Circuit state v(t), at t = 0, ∆t, 2∆t, . . . ≤ T . 1: t ← 0 2: v(t) ← v 0 3: repeat 4: t ← t + ∆t 5: calculate new v(t) 6: until t ≥ T A. Circuit model
The circuit is represented as a set of nodes and components, N and C, respectively. Each component is connected between two nodes, and each node has at least two components connected on it. The circuit is well defined if there are components connecting all the nodes in a closed chain.

For each non-ground node n h ∈ N , we define a set of influencers I h ⊂ N , that groups all the nodes separated from n h by exactly one component. Some authors call these the fanin nodes [START_REF] Newton | Relaxation-based electrical simulation[END_REF].

Consider the circuit in Fig. 2 The iteration is performed in two steps. In the first step, every component reads, from each of its terminal nodes, the voltage v and the sum of current injections φ. Using these values it calculates a new voltage at each node, ψ, and a new current injected to each node, i. Finally it returns these values to the respective nodes. Figure 3 illustrates this procedure for c 1 in our example, which is connected between n B and n C . In the second step every node reads, from each one of the components connected to it, the voltage ψ and the current injection i. Using these values it computes the actual voltage v and the sum of current injections φ. Finally it returns these values to all those components to restart the cycle. Figure 4 illustrates this procedure for n C in our example, whose set of influencers is

I C = {n A , n B , n D }.

B. Component and node models

Now we will enter the boxes of components and nodes and describe how their inputs are converted into outputs.

1) Component:

We represent all components as a voltage source in series with a resistance. This universal component model allows us to describe the behaviour of any of the following: voltage source, resistor, inductor and capacitor. The benefits of using a universal component model will be highlighted later in Section IV, when we discuss implementation issues.

Consider c hk , a component connected between n h and n k . For this component, we define two parameters: vhk , corresponding to the voltage source, measured from n h to n k , and r hk , corresponding to the resistance. Note that vhk = -v kh and r hk = r kh . Next we show how the parameters vhk and r hk are obtained depending on the specific component type.

a) Voltage source: This component consists of an independent voltage source β of maximum current µ. Here, vhk takes the value β, and the resistance r hk is set to a value sufficiently small to ensure that the voltage drop across its two terminals is negligible, i.e., lower than a positive, small real number ǫ:

0 ≤ r hk µ < ǫ. ( 13
)
We achieve this by setting r hk = ǫ 2µ . b) Resistor: This component consists of a constant resistance ̺. Here vhk is set to zero and r hk takes the value ̺.

c) Inductor: This component consists of a constant inductance λ. From Dommel [START_REF] Dommel | Digital computer solution of electromagnetic transients in single and multiphase[END_REF], we know that the dynamic behaviour of an inductance during a time interval [t, t + ∆t], can be emulated by an equivalent circuit of a voltage source in series with a resistance. This gives the values of r hk and vhk in this case:

r hk = 2λ ∆t , (14a)
vhk (t) = v k (t) -v h (t) -r hk i hk (t), ( 14b
)
where t is the simulation clock. d) Capacitor: A constant capacitance γ. We follow the same procedure as for the inductor. Now:

r hk = ∆t 2γ , (15a)
vhk (t) = v h (t) -v k (t) + r hk i hk (t), ( 15b
)
where t is the simulation clock.

v k (t) -v h (t) -r hk i hk (t) Capacitor ∆t 2λ v h (t) -v k (t) + r hk i hk (t)
Table II synthesizes the expressions of vhk and r hk for each component type. Once these parameters are obtained, the outputs of c hk are computed as follows:

ψ hk = v h + r hk φ h , (16a)
i hk = v k -v h -vhk r hk . ( 16b
)
These equations correspond to the Ohm's Law and the Kirchhoff's Current Law applied on c hk , n h and n k . (See Appendix B for the deduction of these equations.) Note that when the sum of current injections on n h is equal to zero, i.e., φ h = 0, we have from equation (16a) ψ hk = v h . In other words, the voltage according to the influencer k is the same as the voltage according to the node itself. This corresponds to a condition of convergence at component-level.

2) Node: Consider now n h and its set of influencers I h . For this node, we define positive real numbers z h and w hk , k ∈ I h such that

z h + k∈I h w hk = 1. (17)
These numbers correspond to the weights of a convex combination that will be used to update v h . In Subsection III-C we show that there is a condition to respect when selecting these weights, in order to ensure convergence.

Once the weights are chosen, the outputs of n h are determined as follows:

v ′ h = z h v h + k∈I h w hk ψ hk , (18a)
φ h = k∈I h i hk . (18b)
Equation (18a) states that a new value of the voltage is computed as a convex combination of its own value, available from the previous iteration, and the values according to all its influencers. Equation (18b) is self-explanatory.

Note that when all the influencers agree on the value of the voltage at n h , i.e., ψ hk = v h , ∀k ∈ I h , we have from equation (18a) v ′ h = v h . This corresponds to a condition of convergence at node-level. When this condition is satisfied at all the nodes, the iteration converges at system-level. In the next subsection we will determine whether this situation is ever reached, and which assumptions need to be made in order to ensure convergence, regarding the theoretical background presented earlier in Section II.

C. Convergence study

Let n G be the ground-node and N ′ the set of all nodes in the circuit minus the ground, i.e., N ′ = N \ {n G }. Let m = |N ′ |, and v ∈ R m be the vector of voltages at all non-ground nodes. We define the function f : R m → R m as follows:

f h (v) = v h + k∈I h j∈I h w hk r hk v j -v h -vhj r hj , h ∈ N ′ . ( 19
)
This expression is obtained by combining equations (16a), (16b), ( 17), (18a) and (18b) and defining v ′ = f (v). Then, f (•) becomes the iteration function. This is a linear application, i.e.,

f (v) = Lv + b, (20)
where L and b are given by:

l hh = 1 - k∈I h j∈I h w hk r hk r hj , l hj =    k∈I h w hk r hk r hj if j ∈ I h , 0 otherwise. b h = - k∈I h j∈I h w hk r hk vhj r hj . (21)
The following theorem establishes that a totally asynchronous iteration associated to f (•) is not certain to converge. Theorem 3. Let v k , k ∈ N be an asynchronous iteration associated to f (•), defined in [START_REF] Dommel | Digital computer solution of electromagnetic transients in single and multiphase[END_REF]. Then, there exist an update function U(•), a delay function d(•) and an initial guess v 0 such that v k does not converge to a fixed point. Proof. To apply Theorem 1 in Section II, we just need to prove that ρ(|L|) ≥ 1, where L is given by [START_REF] Hong | An analytical model for a gpu architecture with memory-level and thread-level parallelism awareness[END_REF]. Now, as the sum of the elements in any row of L is equal to 1, then 1 is an eigenvalue of L and ρ(L) ≥ 1. And by applying Theorem 5 in Appendix A, ρ(|L|) ≥ 1.

This means that to ensure convergence we need further assumptions. The following theorem establishes a sufficient condition for convergence. Theorem 4. Let v k , k ∈ N be a partially asynchronous iteration associated to g(•) defined by:

g(v) = (1 -α)v + α(Lv + b), ( 22
)
where 0 < α < 1 and L and b are given by [START_REF] Hong | An analytical model for a gpu architecture with memory-level and thread-level parallelism awareness[END_REF]. Let the weights z h and w hk , k ∈ I h satisfy 1 -

k∈I h j∈I h w hk r hk r hj ≥ 0, ∀h. ( 23
)
Then v k converges to v * , fixed point of f (•).

Proof. We will use Theorem 2 in Section II. Since we are assuming [START_REF] Anderson | LAPACK Users' Guide[END_REF], we have l hh ≥ 0 (note that l hk ≥ 0 by definition) and then n k=1

|l hk | = n k=1 l hk = 1, ∀h.
So L verifies the assumption in [START_REF] Kapre | Parallelizing sparse matrix solve for spice circuit simulation using fpgas[END_REF]. It remains to prove that L is irreducible. According to Theorem 6 in Appendix A, we need to look at the digraph D(L) and see if it is strongly connected. This digraph can actually be obtained from the circuit diagram, as follows. For each node in the circuit, we place a vertex in D. Next, for each component we place a pair of anti-parallel arcs, connecting the two vertices corresponding to its terminals. Finally, we place an arc from each vertex to itself. Figure 5 illustrates this transformation process. As result, the directed arcs going into each vertex in the graph identify that node's influencers. And as the circuit is well defined, this graph is strongly connected. Fig. 5: Equivalence between the circuit diagram and the digraph of the application matrix L. Some authors call this the dependency graph [START_REF] Newton | Relaxation-based electrical simulation[END_REF].

IV. IMPLEMENTATION AND TESTS

Now that we have described our model, we can complement Algorithm 1 in Section III with the details of the proposed intrinsically parallel iteration. Let C ′ ⊂ C be the set of all inductors and capacitors in the circuit. Algorithm 2 details the implementation of the method using synchronous iterations. The totally asynchronous version can be obtained by eliminating the instruction in line 8; the partially asynchronous, by replacing that instruction by "synchronize processors every d iterations." for all h ∈ C ′ in parallel do update vh (t)

Algorithm
v i ← f i (v)
14:

end for 15: until t ≥ T Note that in instruction 7 of Algorithm 2 all the processors update their coordinate of vector v, by evaluating the iteration function, F (). As we are using a universal component model, the arithmetic operations involved in this function evaluation are actually the same for all processors; only the data differs. This is an ideal situation for the implementation of the model on architectures of the class SIMD (Single Instruction Multiple Data) such as GPU, where such configuration allows for the most data-level parallelism to be exploited.

Next we will present different implementations of Algorithms 1 and 2 along with the results of some tests. We use GPU to implement and run our algorithm, and CPU to implement and run the algorithms from the classic approach. (For a detailed review of GPU architecture, see reference [START_REF] Hong | An analytical model for a gpu architecture with memory-level and thread-level parallelism awareness[END_REF].)

Our computing environment is synthesized in Table III.

A. Implementation alternatives 1) Intrinsically parallel, synchronous (IPS):

We first considered Algorithm 2 and implemented it in CUDA C++ [START_REF] Nvidia | Compute unified device architecture programming guide[END_REF], in order to run it on our Nvidia GPUs. Generally speaking, CUDA is a platform and language for General Purpose GPU programming (GPGPU) allowing to declare and specify kernels, which are executed concurrently by thousands of threads on the GPU. The threads are grouped into blocks and distributed to several CUDA cores, where they can make use of the several GPU resources (e.g., parallel floating-point units, different memory layers, etc.)

The IPS implementation revolves around two classes, Component and Node; both of them are implemented as "structure of arrays," so the program can benefice of coalesced memory accesses (i.e., adjacent GPU threads reading or writing consecutive locations in GPU memory), which is a standard lever for performance. The implementation has two kernels, iterate() and update(), for performing instructions 7 and 13 in Algorithm 2, respectively. These kernels operate at component-level, i.e., the program launches as many threads as components in the circuit and each thread performs the calculations associated to one component. In this way we have a constant number of operations (each component is linked to exactly two nodes), as opposed to a variable number (each node could be linked to any arbitrary number of components) per iteration. This improves the regularity of our solution which is another lever for performance.

Note that CUDA does not provide any explicit block-level mechanism of synchronization. This means that threads from different blocks can get out of phase when running a sequential loop. Therefore in order to ensure synchronization, we let the CPU control the loop in instructions 6 to 9. At every iteration of the loop, the CPU launches the iterate() kernel onto the GPU, so the GPU computes a new approximation of the bus voltages and evaluates the condition of convergence. The CPU waits for the kernel to finish before progressing into the next iteration. The synchronization step in line 8 is implicit in this model.

2) Intrinsically parallel, partially asynchronous (IPPA): Next we considered the modified version of Algorithm 2, in which synchronization between processors is only enforced after a given number of asynchronous iterations. This number is passed as a parameter to the kernel, so the GPU can iterate that many times before returning the control to the CPU. When this parameter is set to one, we obtain the IPS implementation, i.e., the program performs only one (asynchronous) iteration before synchronizing all processors.

3) Classic approach: Last, for fair comparison purposes, we implemented two versions of the classic circuit simulator (see Algorithm 1 in Section III) that relies on matrix factorization to compute the bus voltages. One of these versions utilizes the dense solver LAPACK [START_REF] Anderson | LAPACK Users' Guide[END_REF]; the other one uses the optimized sparse solver KLU [START_REF] Davis | Algorithm 907: Klu, a direct sparse solver for circuit simulation problems[END_REF]. In both cases, the matrix is fully re-factorized at each time step to simulate the case in which there are nonlinearities in the component models.

Both implementations run on single-core CPU. However, it is worth noting that some parallel versions of LAPACK optimized for GPU can yield a speed-up of one order of magnitude [START_REF] Tomov | Dense linear algebra solvers for multicore with gpu accelerators[END_REF]. In KLU, in turn, parallelization seems constrained by the Gilbert-Peierls phase of the algorithm [START_REF] Gilbert | Sparse partial pivoting in time proportional to arithmetic operations[END_REF].

B. Tests and results

Next, in order to evaluate our method (IPS and IPPA implementations) in comparison to the classic approach to circuit simulation (LAPACK and KLU implementations), we randomly generated test circuits with sizes from 500 to above 10,000 nodes. Then we asked our programs to simulate one second of real time operation of these circuits, with a time step of 5 milliseconds, and collected several performance figures. (By default, we use the GTX 680 GPU.)

To generate the test circuits we developed a dedicated C++ program, which uses the Boost library implementation of the Mersenne twister [START_REF] Siek | Boost random number library[END_REF] as a pseudo-random number generator. The strategy used in generating the circuits is the following: first, we generate a random circuit of 16 components that we call a cluster; second, we increase the number of nodes gradually by copying this cluster several times and linking random nodes between the copies; third and final, a random node in every cluster is linked to ground. The idea behind this method is to really isolate the effect of the problem size from all the other factors that could impact on performance (e.g., topology, numerical values, etc.) A similar approach for measuring the performance of the Spice simulator is presented in [START_REF] Newton | Relaxation-based electrical simulation[END_REF].

1) IPS v/s IPPA: Our first experiment consisted in measuring the effect of asynchronicity over performance, by comparing the IPS and IPPA implementations.

Figure 6(a) shows the execution time of IPS and IPPA on four of the test circuits. On the horizontal axis, we have the number of asynchronous iterations performed by the kernel in IPPA (the first point in these curves corresponds to the output of IPS, i.e., the case where the number of asynchronous iterations is set to one). The different curves represent different values of the problem size, n. Clearly, there is a benefit in allowing asynchronism. Furthermore, in this experiment all the four test circuits behave in a very similar way: the simulation time decreases as the number of asynchronous iterations starts increasing, until a certain point after 20 iterations where the gains seem to stall.

Figure 6(b) shows, for the same experiment as before, the average number of iterations executed per time step. Here we observe that the more asynchronous iterations the model is allowed to perform, the more total iterations are needed to converge. In other words, asynchronism has a negative impact on the rate of convergence, which is consistent with the theory exposed in Section II. However, even if that is the case, the overall effect of asynchronism is positive as the cost of these 'extra' iterations is compensated by having to perform fewer global synchronization steps on the CPU side. 2) IPPA v/s KLU and LAPACK: Next we compared IPPA with the classic circuit simulator, represented by the LAPACK and KLU implementations. To maximize performance based in our previous observations, we set the number of asynchronous iterations in IPPA to 40 (see Figure 6(a)).

Figure 7 shows the execution time on a logarithmic scale of IPPA, LAPACK and KLU as a function of the problem size, n. We observe that the IPPA implementation is faster than both the LAPACK and KLU implementations on almost all the considered range. However, as n grows, the simulation time increases faster in IPPA than in KLU. This last situation is most certainly an effect of our implementation not being tuned to yield maximum GPU performance just yet. To illustrate how much room is there for improvement, Fig. 8 presents the IPPA simulation time (left vertical axis) along with the average number of iterations performed per time step (right vertical axis). The horizontal axis contains the problem size, n. Whereas the simulation time grows exponentially with the problem size, the number of iterations seems to be unaffected by n. In other words, it is only the concurrent access to shared resources, which become scarce as the occupancy grows, that prevents maximum performance to be achieved and increases the simulation time more than it can be explained by the amount of work effectively performed by each thread. If the application was tuned to keep every thread working at all times, then it would be only a matter of having enough computing cores for the simulation time to become independent from the problem size. Of course this may be impossible to achieve in practice, as the resources of any computing architecture are essentially limited. However, there are many techniques that can be implemented in order to 'mask' the scarcity of resources as the occupancy grows [START_REF] Nvidia | Compute unified device architecture programming guide[END_REF].

As it is implemented today, we can already observe how the performance of our application scales with the amount of available resources. Figure 9 shows the simulation time of IPPA on different GPU architectures. Notice how, as the number of cores in the architecture grows, the time needed to complete the analysis decreases accordingly. Additionally, the speedup is higher for higher problem sizes.

Our current implementation of IPPA uses GPU's global memory, which is off-chip, to store and load all the relevant data during the course of the iteration. The advantage of doing so is that all the threads can access the data at all times. However, GPU also provides on-chip memory which is only shared by threads within the same thread-block; the access to this shared memory is about one hundred times faster than the access to global memory. Thus, we can conceive an optimization strategy where we assign subsets of stronglyconnected components in the circuit, to threads within the same block on the GPU. Then we let these threads perform a number of local iterations on shared memory, before committing their results to global memory for their integration with those from other blocks. This technique will reduce the number of accesses to global memory by a factor equal to the number of iterations performed locally; however, it might also damage the convergence of our iteration. The full study and implementation of this idea will be addressed in future works.

V. CONCLUSION

We have presented a model for time-domain transient circuit simulation which successfully competes in terms of performance with available state-of-the-art solutions. Furthermore, our approach is intrinsically parallel so it can be efficiently implemented on today's parallel machines. Our tests using a GPU implementation of the model show that the amount of parallel iterations needed to compute the result is independent of the number of nodes in the circuit. However, in order to fully exploit this interesting property, the current implementation must be tuned to achieve maximum utilization of the GPU architecture.

The optimization of the current GPU implementation according to the previous point, as well as the implementation of the proposed method on other parallel architectures, will be addressed in future works. We will also expand our model and investigate convergence issues and applicability of the proposed solution in the non-linear case.

APPENDIX A MATHEMATICAL REMINDERS

In this section we present some classical definitions and mathematical results, for the sake of completion.

Definition 4 (Spectral radius). Let A ∈ C n×n with eigenvalues λ i , 1 ≤ i ≤ n. Then, ρ(A) := max

1≤i≤n |λ i |, (24)
is the spectral radius of A.

Theorem 5 (Bound for the spectral radius [START_REF] Horn | Matrix analysis[END_REF]). Let A ∈ C n×n , and let |A| be the matrix of absolute values of the entries of A. Then ρ(A) ≤ ρ(|A|).

Definition 5 (Reducible matrix). Let A ∈ C n×n . A is called reducible if there exists a permutation matrix P such that

P T AP is of the form     A 1 A 12 0 A 2     ,
where A 1 and A 2 are square matrices of size at least 1. If A is not reducible, then A is called irreducible.

Definition 6 (Digraph of a matrix). Let A = (a ij ) ∈ C n×n . The digraph D with vertices 1, . . . , n, in which there is an arc (i, j) if and only if a ij = 0, is called the digraph of A. We note D = D(A).

Definition 7 (Strong connection). Let D be a digraph and i, j any pair of vertices. If there are directed walks from i to j and from j to i, then D is strongly connected.

Theorem 6 (Interpretation of irreducibility in terms of the digraph of the matrix [START_REF] Horn | Matrix analysis[END_REF]). Let A ∈ C n×n . Then A is irreducible if and only if D(A) is strongly connected.

APPENDIX B DEDUCTION OF THE COMPONENT MODEL EQUATIONS

Assuming the steady-state is attained, in which case v h = ψ hk , we apply Ohm's Law on c hk :

ψ hk = v h = v k -r hk i hk -vhk . ( 25
)
Next we apply Kirchhoff's Current law on n h :

φ h = j∈I h i hj = j∈I h j =k i hj + i hk = 0 (26)
We clear i hk from [START_REF] Siek | Boost random number library[END_REF] and replace it in [START_REF] Gilbert | Sparse partial pivoting in time proportional to arithmetic operations[END_REF]:

ψ hk = v k + r hk j∈I h j =k i hj -vhk . ( 27
)
Expanding, and applying [START_REF] Gilbert | Sparse partial pivoting in time proportional to arithmetic operations[END_REF]:

ψ hk = v k + r hk j∈I h i hj -r hk i hk -vhk = v h + r hk j∈I h i hj . (28)
Finally, taking φ h from [START_REF] Siek | Boost random number library[END_REF] and replacing it in (28) we obtain the first formula:

ψ hk = v h + r hk φ h . (29)
And clearing i hk from [START_REF] Gilbert | Sparse partial pivoting in time proportional to arithmetic operations[END_REF] we obtain the second formula:

i hk = v k -v h -vhk r hk . (30)

9 Fig. 1 :

 91 Fig. 1: Bound on the convergence rate as a function of the problem size, for different values of the parameter σ.

Fig. 2 :

 2 Fig. 2: Example circuit.

Fig. 3 :

 3 Fig. 3: Component-level procedure.

Fig. 4 :

 4 Fig.4: Node-level procedure.

  Digraph of matrix L.

Fig. 6 :

 6 Fig. 6: Effect of asynchronism on performance: (a) simulation time and (b) average number of iterations per time step, as functions of the number of asynchronous iterations allowed.

Fig. 7 :

 7 Fig. 7: Related performance of IPPA, LAPACK and KLU.

Fig. 8 :Fig. 9 :

 89 Fig. 8: Simulation time and number of iterations performed by IPPA.

TABLE I :

 I List of variables and parameters of the model.

	Symbol	Type	Definition
	N	Fixed	Set of circuit's nodes.
	C	Fixed	Set of circuit's components.
	I h	Fixed	Set of influencers of node h.
	v h	Variable	Voltage at node h.
	φ h	Variable	Sum of current injections on node h.
	ψ hk	Variable	Voltage at node h according to the influencer
			k.
	i hk	Variable	Current flowing from node k to node h.
	r hk	Fixed	Resistance of the component situated between
			nodes h and k.
	vhk	Fixed	Voltage source of the component situated be-
			tween nodes h and k.
	z

h Fixed Weight within node h of its own voltage. w hk Fixed Weight within node h of the voltage determined by the influencer k.

TABLE II :

 II Model parameters per component type.

	Component	r hk	vhk
	Voltage source	ǫ 2µ	β
	Resistor	̺	0
	Inductor	2λ ∆t

TABLE III :

 III Computing environment used in the tests.

	Hardware	Clock rate (GHz)	Number of cores
	Xeon X560 CPU		2.67	12 (1 used)
	GeForce GTX 560 GPU		1.62	336
	GeForce GTX 480 GPU		1.40	480
	GeForce GTX 680 GPU		1.06	1,536
	Software	Version
	GCC		4.8.2
	CUDA	6.0
	LAPACK	3.5.0
	KLU		1.2.1